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I. 


INTRODUCTION 


A.  PURPOSE 

The  purpose  of  this  thesis  was  to  study  the  effects  of  various  geometries  on  the 
reflected  shock  conditions  for  producing  high  local  temperature  and  pressure  regions.  The 
nominal  geometry  was  representative  of  that  existing  in  a  valve-less  pulse  detonation 
engine  (PDE)  developed  by  the  Naval  Postgraduate  School  (NPS).  Currently,  the  NPS 
PDE  uses  a  small  auxiliary  combustor  known  as  an  initiator,  which  operates  on  an 
oxygen-enriched  fuel/air  mixture  to  rapidly  and  reliably  generate  a  detonation  wave, 
which  then  re-initiates  the  less  sensitive  fuel/air  mixture  in  the  combustor.  Since  the 
addition  of  oxygen  is  treated  as  “fuel”  for  specific  impulse  (ISP)  calculation  and  reduces 
the  overall  system  performance,  it  is  undesirable  on  any  propulsion  platform  in  terms  of 
performance  and  practically. 

The  benefits  of  modified  combustor  geometry  are  to  eliminate  the  requirement  for 
auxiliary  oxygen  in  the  initiator  and  only  utilize  an  exiting  fiiel/air  detonation  wave  from 
the  initiator  into  the  main  combustor.  The  modified  geometry  should  result  in  the 
generation  of  multiple  regions  of  high  temperature  due  to  the  reflected/focused  shock 
conditions  from  its  geometry,  and  they  subsequently  should  allow  the  initiation  of  a 
detonation  of  the  fiiel/air  mixture  in  the  main  combustor. 


B.  INTRODUCTION  TO  PULSE  DETONATION  ENGINE 

The  interest  in  PDEs  has  increased  dramatically  in  recent  years  due  to  their  high 
theoretical  performance  and  wide  range  of  potential  applications.  PDEs  use  detonation 
waves  that  propagate  through  a  premixed  fiiel/air  mixture  to  produce  large  increases  in 
chamber  pressures  and  temperatures.  The  expansion  of  these  products  produces  thrust. 
They  offer  the  potential  for  high  performance  because  the  rapid  detonation  process  results 
in  a  very  efficient  combustion  process  with  high  operating  frequencies  and  combustor 
chamber  pressures.  They  offer  the  potential  of  good  thrust  characteristics  across  subsonic 
to  supersonic  flight  regimes.  They  have  the  potential  to  operate  at  very  high  energy 


1 


densities,  allowing  the  use  of  simple  and  compact  combustor  designs.  They  also  offer 
instant  throttleability  due  to  cycle-to-cycle  control. 

The  valve-less  PDE  developed  at  the  NPS  uses  a  continuous  air  flow  design, 
which  eliminates  the  need  for  valves  to  supply  fuel/air  mixture  to  the  combustor.  The 
initiator  operates  on  an  oxygen-enriched  fuel/air  mixture  to  rapidly  and  reliably  generate 
a  detonation  wave,  which  re-initiates  the  less  sensitive  fuel/air  mixture  in  the  combustor. 
Without  a  large  valve  on  the  air  flow,  a  convenient  flow  path  is  allowed  to  rapidly  fill, 
detonate,  and  purge  the  combustion  chamber  at  rates  up  to  100  Hz  cycle  frequency. 

C.  INTRODUCTION  TO  COMPUTATIONAL  FLUID  DYNAMIC  TOOLS 

Computational  Fluid  Dynamics  (CFD)  is  a  powerful  tool  for  fluid  flow  analysis. 
The  CFD  process  involves  many  primary  and  secondary  programs  used  on  Silicon 
Graphic  Iris  (SGI)  workstations  at  the  NPS.  The  lists  of  programs  used  for  this  thesis  are 
as  follows: 

•  GRIDGEN  -  Grid  generation  program 

•  GRIDED  -  Grid  editing  program 

•  OVERFLOW  -  Flow  solver  program 

•  FAST  -  Flow  solution  visualization  program 

The  interpretation  of  CFD  results  requires  much  scientific  judgment  and 
knowledge.  The  CFD  results  can  be  improved  by  modifying  the  grid  mesh,  boundary 
conditions,  numerical  differencing  scheme,  time  step  and  more. 

The  OVERFLOW  flow  solver  program  is  based  on  prefect  gas  law  algorithm; 
solution  will  be  stable  for  flow  less  than  Mach  5  and  it  is  not  capable  of  handling 
chemical  reaction.  Due  to  its  limitations,  the  model  is  not  capable  of  simulating  the 
operating  conditions  of  the  actual  PDE.  But  it  is  believed  to  be  adequate  to  study  the 
efforts  of  combustor  geometry  on  shock  reflecting/focusing  to  generate  regions  of  high 
vorticity  and  subsequently  regions  of  high  temperature  also  know  as  “hot  spots”, 
transmitted  from  the  initiator  into  the  main  combustor  through  a  weaker  supersonic  shock 
wave. 
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II.  DETONATION  PHYSICS 


A.  INTRODUCTION 

This  chapter  discusses  the  physics,  properties,  structure  and  development  of 
detonation.  Detonation  is  a  supersonic  combustion  wave  that  typically  propagates  at  a 
few  thousand  meters  per  second,  closely  coupled  to  a  thin  flame  front.  Detonation  waves 
travel  in  this  manner,  which  are  locally  three-dimensional  and  non-steady.  Short 
deflagration  to  detonation  transition  (DDT)  length  is  desirable  for  propulsion  system 
design  and  can  be  achieved  by  increasing  the  chemical  reaction  rate  of  premixed  gases 
through  turbulence  and  increasing  the  temperature  through  shock  reflection. 


B.  EXPLOSION,  DFLAGRATION  AND  DETONATION 

The  chemical  reactions  relating  to  the  combustion  of  premixed  gases  generally 
fall  into  one  of  the  three  categories  namely  explosion,  deflagration  and  detonation.  The 
three  categories  can  be  further  described  as  follows.  The  chemical-reaction  zone  is  often 
called  the  “flame  zone”,  “flame  front”  or  “reaction  wave”,  which  often  emits  visible  light 
[Ref  1].  There  are  generally  two  types  of  flames  namely  premixed  flames,  where 
reactants  are  perfectly  mixed  before  chemical  reaction  and  diffusion  flames,  which 
reactants  diffuse  and  mix  into  each  other  during  the  chemical  reaction  [Ref  1]. 

1.  Explosion 

Explosion  is  a  term  that  corresponds  to  a  rapid  heat  release  or  pressure  rise  but  it 
does  not  require  the  passage  of  a  combustion  wave  through  the  exploding  medium. 
Certain  fuel  and  oxidizer  gas  mixtures  will  not  propagate  a  burning  zone  or  combustion 
wave  and  are  outside  the  flammability  limits  of  the  explosive  gas  [Ref  2]. 

2.  Deflagration 

Deflagration  is  a  combustion  wave  propagating  at  subsonic  speed  and  is  the  most 

common  type  of  combustion  [Ref  1].  Deflagration  flame  speeds  at  the  order  of  one  or 

more  meters  per  second  for  most  typical  hydrocarbon/air  mixtures.  The  propagation  of  a 

deflagration  flame  is  governed  by  laminar  or  turbulent  diffusion  of  unbumed  gases  ahead 
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of  a  flame  and  burned  gases  behind  the  flame.  Deflagrations  produce  small  decreases  in 
pressure  and  can  be  modeled  as  nearly  isobaric  or  constant  pressure  processes  [Ref.  1]. 
Some  examples  of  deflagrations  are  gas  stove  flames  and  a  camp  Are  flame. 


3.  Detonation 

Detonation  is  a  supersonic  combustion  wave  that  typically  propagates  at  a  few 
thousand  meters  per  second  relative  to  an  unbumed  fuel/air  mixture  [Ref  1].  Detonations 
are  much  more  energetic  than  deflagrations  and  produces  large  overpressures  because  the 
combustion  region  is  closely  coupled  to  a  supersonic  shock  wave  [Ref  3]. 

A  detonation  wave  compresses  a  fluid,  increasing  its  pressure,  density  and 
temperature.  The  rapid  compression  and  heating  of  the  reactants  greatly  accelerates  the 
combustion  process.  Detonations  can  be  approximated  as  normal  shock  waves  that 
initiate  and  are  closely  coupled  to  a  thin  flame  front  or  combustion  region.  A  one¬ 
dimensional  model  of  such  a  structure  is  known  as  a  Zel’dovich,  von  Neumann  and 
Doring  (ZND)  model,  which  is  further  described  later  in  this  chapter.  Due  to  the  high¬ 
speed  nature  of  a  detonation  wave,  a  detonation  closely  approximates  a  constant  volume 
combustion  process,  thereby  producing  a  lower  entropy  rise  of  the  working  fluid  versus  a 
deflagration  combustion  mode.  The  result  is  an  improved  thermodynamic  efficiency  [Ref 
1]- 


4.  Qualitative  Differences  between  Detonation  and  Deflagration 

A  schematic  diagram  of  a  one-dimensional  combustion  wave  is  shown  in  Figure  1 
[Ref  1].  In  this  reference  frame,  the  unbumed  reactants  move  towards  the  stationary 
combustion  wave  at  velocity,  ui,  and  the  burned  reactants  move  away  from  the 
combustion  wave  at  velocity,  U2.  The  difference  between  these  two  velocities  can  be  large 
depending  on  whether  the  combustion  process  occurs  as  deflagration  or  detonation. 
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stationary  Combustion  Wave 


Figure  1 .  Schematic  of  Stationary  One-Dimensional  Combustion  Wave. 


A  comparison  of  thermodynamic  properties  between  deflagration  and  detonation 
is  shown  in  Table  1  [Ref.  1].  The  terms  u,  c,  p,  T  and  p  are  velocity,  sonic  velocity, 
pressure,  temperature  and  density  respectively. 


Ratio 

Detonation 

Deflagration 

Ui/Ci 

5-10 

0.0001-0.03 

U2/U1 

0.4-0. 7  (deceleration) 

4-6  (acceleration) 

P2/P1 

13-55  (compression) 

==0.98  (slight  expansion) 

T2/T1 

8-21  (heat  addition) 

4-16  (heat  addition) 

P2/P1 

1. 7-2.6 

0.06-0.25 

Table  1 .  Qualitative  Differences  between  Detonation  and  Deflagration  in  Gases. 


C.  HUGONIOT  CURVE  AND  CHAPMAN-JOUGUET  POINTS 

A  detonation  wave  can  be  approximated  as  a  normal  shock  wave  that  is  closely 
followed  by  a  deflagration.  However,  detonation  waves  are  actually  complex,  oscillatory 
phenomena  with  three-dimensional  cellular  structure  [Ref  3]. 

Detonation  waves  can  be  modeled  as  discontinuities  in  the  flow  of  an  ideal  gas  at 
which  heat  addition  occurs  [Ref.  4].  This  treatment  of  detonation  waves  is  a  reasonable 
approximation  for  studying  engine  performance  since  the  significant  changes  in  the  fluid 
properties  occur  over  a  very  short  distances.  The  analysis  is  based  on  the  conservation 
equations  for  mass,  momentum,  and  energy  across  the  detonation  in  a  constant  cross 
section  tube  and  on  the  equation  of  state  [Ref.  3]. 

Conservation  of  Mass: 


yOjMj  =  =  m 


(1) 
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Conservation  of  Momentum: 


+  P]U^  — 


(2) 


Conservation  of  Energy: 


+  2  ^  “  ^^^2  +  2  ^2^ 

/  1  2  _,  1  2 

/Zi  +  —  t/i  +  ^  —  /Z2  +  “t/2 
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^  r; 


vr-ly 


R 


where 


(3) 

(4) 


subscript  1  refers  to  state  of  unbumed  gas  upstream  of  the  detonation  wave 
subscript  2  refers  to  state  of  burned  gas  immediately  downstream  of  detonation 


p  =  density,  v  =  specific  volume  =  1/p 
u  =  one-dimensional  velocity 
h  =  enthalpy  of  the  gas  (per  unit  mass) 
q  =  heat  addition  (per  unit  mass) 

Combining  the  conservation  of  mass  (1)  and  conservation  of  momentum  (2) 
equations  produces  the  Rayleigh-line  relation  [Ref  1]  as  defined  in  equations  (4)  and  (5)  in 
terms  of  Mach  number.  Since  (piui)^and  (p2U2)^  are  always  positive,  the  terms  (P2-P1)  and  (V2-V1) 
must  be  of  opposite  sign. 


P,-P,=  p,u^  -  p^u^ 
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(5) 


A  deflagration  produces  a  net  decrease  in  both  pressure  and  density  as  a  result  of 
combustion  process.  Therefore  deflagrations  are  classified  as  being  weak  because  they  are  nearly 
isobaric  and  result  in  only  small  decrease  in  pressure.  However,  a  detonation  produces  a  net 
increase  in  both  pressure  and  density  because  the  shock  wave,  which  precedes  the  flame  region  of 
the  detonation  wave,  both  compresses  and  pressurizes  the  un-reacted  gas  prior  to  the  heat  addition 
process  [Ref  3]. 

Combining  and  rearranging  die  conservation  of  mass  (1),  conservation  of  momentum  (2) 
and  conservation  of  energy  (3)  equations  forms  the  Hugoniot  equations  for  flow  in  a  constant  cross- 
sectional  tube  with  heat  release.  Equations  (6)  and  (7)  relate  the  thermodynamic  properties  upstream 
and  downstream  of  a  combustion  region  [Ref  1  ] . 

K-h^+q  =  ^{P,-P^){v^-v,)  (6) 

^i-e2+^  =  ^(^+^2)(n-^2)  (7) 

where 

e  =  internal  energy  of  gas  per  unit  volume 

h  =  enthalpy  =  e+p/p  =  e+pv 

Champan  (1899)  and  Jouguet  (1905)  [Ref  3]  independently  determined  that  the  tangency 
point  of  the  Rayleigh  line  and  the  Hugoniot  curve  for  the  detonation  reaction  products  represented 
the  stable  end  state  for  self-sustaining  detonation  waves.  The  schematic  diagram  depicting  the 
various  sections  of  the  Hugoniot  curve  is  shown  in  Figure  2  [Ref  4].  The  Chapman- Jouguet  (CJ) 
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condition  for  detonation  is  indicated  by  point  B  (P+,  v+).  A  Rayleigh  line  through  the  point  (1,1) 
win  not  intersect  the  Hugoniot  curve  if  the  magnitude  of  its  slope  is  less  than  die  tangent  line’s 
slope,  which  implies  a  minimum  wave  speed,  p,  for  detonations.  There  are  two  possible  detonation 
end  states  at  higher  wave  speeds:  one  on  curve  A-B  representing  strong  detonation  and  another  on 
curve  B-C  representing  weak  detonations.  The  summary  of  the  properties  of  Hugoniot  curves  is 
shown  in  Table  2  [Ref  4]. 

P 


Section  in 
Figure  2 

Pressure 

Ratio 

P=(Poo/po) 

Velocity  and 
density  ratios 

\^(Voo/Vo)=(po/poo) 

Propagation 
Mach  number 

Mo=(vo/af,o) 

Downstream 
Mach  Number 

Moo=(Voo/ae,oo) 

Remarks 

Strong 

detonations 

Line  A-B 

P+<P<00 

niin<v<  v+ 
(vniin>0) 

Mo+^Mq^oo 

Moo<l 

Seldom  observed; 
required  special 
experimental  arrangement 

Upper 

Chapman- 
Jouguet  point 

Point  B 

p=p+ 

(p^l) 

v=v+ 

(V+<1) 

Mo=Mo+ 

(Mo+>l) 

Moo=l 

Usually  observed  for 
waves  propagating  in 
tubes 

Weak 

detonations 

Line  B-C 

Pl<p<p+ 

(P.>1) 

V+<V<1 

Mo+^Mq^oo 

Moo>l 

Seldom  observed;  requires 
very  special  gas  mixtures. 

Weak 

deflagrations 

Line  D-E 

p.<p<l 

V  1<V<  V. 

(v  1>1) 

O^Mq^Mq- 

Moo<l 

Often  observed;  p~l  in 
most  experiments. 

Lower 

Chapman- 
Jouguet  point 

Point  E 

p=p- 

(P-<1) 

V  =  V. 

(V>1) 

Mo=Mo- 

(Mo-<l) 

Moo=l 

Not  observed. 

Strong 

deflagrations 

Line  E-F 

0<p<p. 

V  <V<Vmax 

(v  max<00) 

Momin<Mo<Mo- 

(Momin^O) 

Moo>l 

Not  observed;  forbidden 
by  considerations  of  wave 
structure. 

Table  2.  Properties  of  Hugoniot  Curves. 
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D.  DETONATION  WAVE  PROPERTIES 

A  general  pressure  distribution  of  a  detonation  initiated  at  the  closed  end  of  a  tube 
and  propagating  left  to  right  toward  the  open  end  is  shown  in  Figure  3  [Ref  1]. 


vEjnNeomann 

© 


Unbumed 


I  P2,T  Ps,Ts  Mbctuire 

Figure  3.  Detonation  Chamber  Pressure  Profile. 


The  transitory  von  Neumann  pressure  increase  across  the  leading  shock  of  the 
leading  shock  of  the  detonation  wave  can  be  defined  from  the  conservation  equations  (1), 
(2)  and  (3)  and  the  perfect  gas  relation  as: 

^  =  1  +  ^^(m,^-1)  (8) 

where 

yi  =  specific  heat  ratio  in  unbumed  gas 

Ml  =  Mach  number  of  detonation  wave  relative  to  unbumed  gas 

The  following  relations  can  be  expressed  for  pressure  and  temperature  ratios 
across  a  Chapman- Jouguet  detonation  wave  in  an  ideal  gas: 

^2  _  1  +  fi-^i 

^  l  +  f2 


t;  Wj  ri  ^1 A  1  +  72  y 


where 
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y  =  ratio  of  specific  heats 
m  =  molecular  weight  of  gaseous  mixtures 

Relating  the  speed  of  the  unbumed  gases  just  behind  the  detonation  to  the  speed 
of  the  wave  and  accounting  for  the  initial  velocity  due  to  fuel/air  injection  yields: 


I  72  T2  y 


-1 


(11) 


where 


M2C  =  Mach  number  of  burned  gas  w.r.t.  detonation  chamber  wall 
Mo  =  Mach  number  of  unbumed  gas  w.r.t.  detonation  chamber  wall 
Ml  =  Mach  number  of  detonation  wave  w.r.t.  unbumed  gas 

Mach  number  of  the  burned  gas  relative  to  the  detonation  chamber  wall,  M2C  is  in 
the  high  subsonic  range  for  detonations  initiated  within  tubes,  which  contain  a  fuel/air 
mixture  with  no  initial  velocity.  M2C  can  be  more  than  one  if  the  velocity  of  the  unbumed 
gas  is  not  zero. 

The  von  Neumann  spike  pressure,  Ps,  is  typically  a  large  value.  Due  to  the 
boundary  condition,  the  static  pressure  at  the  closed  end  is  reduced  to  P3  by  a  series  of 
rarefaction  waves  that  forms  in  order  to  ensure  the  normal  velocity  is  zero  at  the  closed 
end.  The  pressure  just  behind  the  detonation  wave,  P2,  is  therefore  reduced  to  the  plateau 
pressure,  P3,  at  the  wall.  P3  can  be  deduced  from  the  following  pressure  ratio: 


^2-1 

2 


M2C 


\ 


2/2 

fl-l 


J 


(12) 


E.  DETONATION  WAVE  STRUCTURE 

1.  ZND  One-Dimensional  Detonation  Wave  Structure 

ZeTdovich  (1940),  von  Neumann  (1942)  and  Doring  (1943)  independently 
assumed  that  the  flow  is  strictly  one-dimensional  and  steady  relative  to  the  detonation 
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front  [Ref.  1],  The  variation  of  physical  properties  through  one-dimensional  ZND 
detonation  wave  is  shown  in  Figure  4  [Ref  1]. 


Shock 


X 


Figure  4.  Variation  of  Physical  Properties  through  a  ZND  Detonation  Wave. 

The  shock  wave  moves  at  the  detonation  velocity,  with  chemical  reaction 
occurring  behind  it  in  a  region  much  thicker  than  a  typical  shock  wave.  It  initially  heats 
the  reactants  to  a  temperature  at  which  they  can  react  at  a  high  enough  rate  for 
deflagration  to  locally  propagate  as  fast  as  the  wave.  Its  thickness  is  usually  very  small,  in 
the  order  of  a  few  mean  free  paths  of  the  gas  molecules.  Thus  the  major  heat  release 
occurs  in  a  substantially  thicker  region  behind  it  [Ref  1].  The  pressure,  temperature  and 
density  at  state  1  increase  sharply  behind  the  shock  at  state  1’.  According  to  Arrhenius 
law  as  shown  below  in  equation  (13),  the  reaction  rate  increases  in  the  region 
immediately  behind  the  shock,  where  the  temperature  is  not  extremely  high.  This  region 
is  called  the  induction  zone  because  the  pressure,  temperature  and  density  profile  are 
relatively  flat  immediately  behind  the  shock  front.  At  the  reaction  zone,  the  gas  properties 
change  sharply  as  the  reaction  rate  increases  drastically  to  high  values  and  reaches 
equilibrium  values  at  the  end  of  the  reaction.  The  distance  between  the  shock  front  and 
the  fully  reacted  state  for  most  fuel/air  mixtures  is  in  the  order  of  1cm  [Ref  1]. 
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The  locus  of  the  reacting  mixture  on  the  Hugoniot  curve  is  shown  in  Figure  4 
[Ref.  1].  Paths  (a),  (b),  (c)  and  (d)  pass  through  the  detonation  wave  from  un-reacted 
state  to  fully  reacted  state,  therefore  they  satisfy  the  conservation  equations.  Path  (a)  is 
unlikely  to  have  sufficiently  fast  reaction  rates  to  sustain  the  wave  whereas  paths  (b)  and 
(c)  represent  mixtures  with  fast  and  slow  chemical  kinetics  respectively.  The  peak 
pressure  behind  the  shock  wave  is  known  as  the  von  Neumann  spike. 


Figure  5.  ZND  Detonation  Structure  on  (p,l/p)  Diagram. 

The  Arrhenius  rate  law  [Ref  1]  is  defined  as: 

(SI 

l,  =  Ae^^  '  (13) 

where 

k  =  rate  of  reaction 
A  =  constant  to  each  chemical 
Ea  =  activation  energy 
Ru  =  universal  gas  constant 


T  =  temperature 
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Therefore  the  chemical  reaction  rate  increases  exponentially  as  the  temperature  of 
the  mixture  increases.  The  constant,  A,  generally  increases  as  P  . 

2.  Multi-Dimensional  Detonation-Wave  Structure 

By  observing  spin  in  limit  mixtures  in  circular  tubes,  they  showed  that  detonation 
waves  might  prefer  to  travel  in  a  manner,  which  is  locally  three-dimensional  and  non¬ 
steady  [Ref  1].  Through  the  smoked- foil  method,  a  propagating  detonation  wave  is  able 
to  write  on  the  walls  due  to  the  presence  of  triple  points  at  the  intersecting  lines  of  three 
shock  waves  as  shown  in  Figure  6  [Ref  1],  A  steady  supersonic  two-dimensional  flow, 
passing  through  a  convergent  ramp  section  from  right  to  left,  produces  a  shock-wave 
pattern,  when  the  ramp  is  very  steep  as  shown  in  Figure  7  [Ref  1].  The  triple  point  is  the 
intersection  point  of  Mach-stem,  incident  and  reflected  shocks.  A  slipstream  is  also 
generated. 


Figure  6.  Smoked-Foil  Record  and  Schematic  Diagram  of  Symmetric  Planer  Interaction. 
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The  Defonttfon  Wove 


Unburned  Region 
(Uniform  Properties) 


Burned  Region 

Uniform  Pressure  a  Velocity 


Bockground  I  Much  Stem  Structure  In  a  Reoctive  Gos 


Figure  7.  Schematic  Diagram  Showing  the  Shock-Wave  Pattern  and  Triple  Point  in  a  Two- 
Dimensional  Supersonic  Flow  Passing  through  a  Convergent  Ramp  Section. 


Cell  length,  Lc,  is  defined  as  the  cell  dimension  along  the  axis  of  detonation 
propagating  and  cell  size,  X,  also  known  as  cell  width,  is  defined  as  the  cell  dimension 
perpendicular  to  the  direction  of  propagation  as  shown  in  Figure  8.  The  chemical 
reactions  are  essentially  completed  within  one  cell  length,  and  X  =0.6  Lc.  The  minimum 
cell  size,  X,  usually  occurs  at  equivalence  ratio,  (p,  of  1 .  Equivalence  ratio,  (p,  is  defined  as 
the  ratio  of  actual  fiiel-oxidant  ratio  to  the  fuel-oxidant  ratio  for  a  stoichiometric  process. 
Cell  size,  X,  is  also  representative  of  the  sensitivity  of  the  mixture.  It  increases  as  the 
sensitivity  of  the  mixture  decreases  [Ref  2].  Figure  9  shows  the  detonation  cell  size 
against  equivalence  ratio  for  common  fuels  [Ref  6]. 


14 


LEGEND 

Fuel  Smoke  Pressure 

Foil  Oscillations 
CH,  ^ 

C4H.,,  • 

CjH,  V  V 

CjH,  □ 

CjH,  A 

H2  0  # 

CjHj  X 

0  1  2  3  4 

EQUIVALENCY  RATIO 

Figure  9.  Cell  Size  Measurements  for  Common  Fuel. 


F.  DEVELOPMENT  OF  DETONATION  WAVE 

Direct  ignition  and  deflagration-to-detonation  transition  (DDT)  process,  also 

known  as  thermal  ignition,  are  primarily  the  two  ways  to  initiate  a  detonation  wave.  DDT 

length,  Xddt,  is  dependent  on  the  type  of  combustible  mixture  used  and  its  equivalence 

ratio,  the  internal  combustor  geometry  and  possible  heat  addition.  Figure  10  shows  a 
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Figure  10.  Deflagration-to-Detonation  Transition  for  Hydrogen/Oxygen  Mixture. 


A  cylindrical  tube  contains  a  combustible  mixture,  which  is  ignited  with  a  spark 
plug  at  the  closed  end  and  the  other  end  opens.  A  combustion  wave  initially  propagates  as 
a  laminar  flame  front  and  rapidly  transit  into  a  turbulent  flame  front,  propagating  at  5  to 
1 0  times  that  of  the  laminar  speed.  This  flame  produces  weak  compression  waves,  which 
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propagate  at  the  local  sonic  velocity  and  increases  the  temperature  of  the  unbumed 
mixture  slightly  after  each  wave  passes.  A  shock  wave  rapidly  develops  as  these 
compression  waves  coalesce.  Based  on  the  local  sonic  velocity  defined  in  equation  (14), 
the  increased  temperature  behind  the  shock  front  increases  the  sonic  velocity  further  [Ref 
5]. 

c  =  4^  (14) 

where 


c  =  local  sonic  velocity 

y  =  ratio  of  specific  heats  for  reactants  and  hot  products  throughout  reaction 
R  =  gas  constant 
T  =  local  temperature 

A  combustion  driven  shock  is  generated  as  the  turbulent  flame  continues  to 
accelerate  through  the  shock-heated  and  compressed  mixture.  The  combustion  wave 
continues  to  produce  compression  waves,  which  propagate  towards  the  shock  and  further 
increase  the  temperature  of  the  unburned  mixture.  When  a  sufficient  temperature  rise 
occurs  between  the  shock  and  the  combustion  front,  a  thermal  explosion  will  occur  in  the 
unbumed  mixture,  producing  strong  secondary  shock  waves  in  all  directions  and  allowing 
the  combustion  zone  to  occur  immediately  after  the  leading  shock  wave  forming  a 
detonation  wave.  The  detonation  will  overdrive  itself  momentarily,  but  then  propagates 
at  about  Mach  5  relative  to  the  unbumed  mixture.  The  reaction  zone  then  becomes 
located  close  behind  the  leading  shock  of  the  detonation  wave.  The  detonation  shock 
wave  preheats  the  unbumed  mixture  to  extremely  temperature  as  it  propagates  through 
causing  the  reactants  to  combust  extremely  fast.  The  equations  governing  the  Mach 
number  and  temperature  ratio  before  and  after  the  shock  in  a  prefect  gas  are  defined  in 
equations  (15)  and  (16)  [Ref  5]. 


M 


2 

2 


M/  + 


(f-1) 


vf-ly 


Ml  -1 


(15) 
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(16) 


2  " 

where 

Ml  =  Mach  no.  of  the  leading  shock 

M2  =  Mach  no.  after  the  shock 

Ti  =  temperature  before  the  shock 

T2  =  temperature  after  the  shock 

Typical,  Xddt  for  most  hydrocarbon/air  combustion  is  ten  tube  diameter  or  more 
for  a  tube  diameter  greater  than  the  cell  size,  X  [Ref  6].  Xddt  must  take  place  within 
about  0.3  m  as  the  practical  lengths  for  propulsion  purposes  are  about  1  m.  Short  Xddt  is 
desirable  for  propulsion  system  design  and  can  be  achieved  by  increasing  the  global 
reaction  rates  through  turbulence  and  increasing  the  temperature  through  shock  reflection 
[Ref  5].  High-speed  turbulent  combustible  mixtures  injected  through  orifices  into  the 
combustor  can  also  increase  the  overall  reaction  rate  as  can  reflected  shock  scenarios  due 
to  combustor  geometry  [Ref  5]. 
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III.  PULSE  DETONATION  ENGINES 


A.  INTRODUCTION 

This  chapter  discusses  the  operating  principles  of  PDEs.  PDEs  have  received 
renewed  interest  during  the  past  10  years.  They  use  detonation  waves  that  propagate 
through  a  premixed  fuel/air  mixture  to  produce  large  chamber  pressures,  which  then 
produce  thrust  [Ref.  3].  They  offer  the  potential  for  high  performance  because  the  rapid 
detonation  process  results  in  nearly  constant  volume  combustion  with  high  operating 
frequencies  and  combustor  chamber  pressures.  They  are  predicted  to  be  very  efficient  and 
offer  good  thrust  characteristics  from  subsonic  to  the  supersonic  flight  regimes.  They 
have  the  potential  to  operate  at  very  high  energy  densities,  allowing  the  use  of  simple  and 
compact  combustor  designs  [Ref  3].  They  also  offer  instant  throttleability  due  to  cycle- 
to-cycle  control. 


B.  STEADY  AND  UNSTEADY  STATE  ENGINES 

Air-breathing  engines  can  be  classified  according  to  the  type  of  combustion 
process  employed  in  the  device.  The  combustion  process  can  be  characterized  as  either 
steady  (quasi-steady)  or  unsteady.  Propulsion  systems  may  be  further  classified  according 
to  whether  a  deflagrative  or  detonative  mode  of  combustion  is  utilized  [Ref  3]. 

The  most  commonly  used  class  of  engines  is  steady  state  combustion  propulsion 
systems.  A  continuous,  constant  pressure,  deflagration  combustion  process,  characterizes 
these  engines.  Examples  of  this  class  of  engines  are  the  turbojets  and  ramjets. 

The  second  class  of  engines  is  the  unsteady  combustion  propulsion  systems. 
Combustion  occurs  in  an  intermittent  manner  and  is  governed  by  the  fuel/air  flame 
speeds  and  combustion  chamber  characteristics.  Unsteady  combustion  propulsion 
systems  can  be  based  on  either  deflagration  combustion  (pulse  jets)  or  detonation 
combustion  (pulse  detonation  engines)  [Ref  3]. 
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C.  ENGINE  EFFICIENCY 

The  constant  volume  detonation  and  constant  pressure  deflagration  processes  can 
be  modeled  as  the  Humphrey  and  Brayton  cycles  respectively.  Figures  11  and  12 
compare  the  temperature,  pressure  and  specific  volume  characteristics  of  the  Brayton  and 
Humphrey  cycles  [Ref  1].  The  Brayton  cycle  (0-1 -4-5-0)  consists  of  two  constant 
pressure  cycle  (1-4  and  5-0)  and  two  isentropic  processes  (0-1  and  4-5).  The  Humphrey 
cycle  is  similar,  except  the  constant  pressure  combustion  process  (1-4)  is  replaced  by  a 
constant  volume  heat  addition  process  (1-2). 


Figure  11.  Pressure- Volume  Cycle  Diagram  (Left). 

Figure  12.  Temperature-Entropy  Cycle  Diagram  (Right). 


The  efficiencies  of  the  constant  pressure  Brayton  and  constant  volume  Humphrey 
cycles  can  be  computed  from  the  pressure-volume  and  temperature-entropy  diagrams  in 
Figures  1 1  and  12.  The  efficiency  of  a  cycle  is  defined  as  the  useful  work  output  divided 
by  the  total  heat  energy  input. 

The  efficiency  of  the  Brayton  cycle  defined  in  equation  (17)  depends  only  on  the 
temperature  change  during  either  of  the  two  isentropic  compression  or  expansion 
processes  (i.e.  Tq/Ti  =  T4/T5)  [Ref  3]. 


^BRAYTON  rr 


(17) 
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The  efficiency  of  the  Humphrey  cycle  [Ref  3]  is  defined  as: 


V HUMPHREY  1  y  rj. 


f  rr.  \ 


T 

^-1 


(18) 


The  efficiency  of  the  Humphrey  cycle  depends  not  only  on  the  isentropic 
compression  temperature  ratio,  Tq/Ti,  but  also  on  the  ratio  of  specific  heats,  y,  and  the 
temperature  change  due  to  the  constant  volume  combustion  (i.e.,  the  detonation 
temperature  ratio  T2/T1). 

The  difference  between  the  Brayton  and  Humphrey  cycle  efficiencies  is  the 
following  Tq/Ti  multiplier: 
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The  value  of  this  expression  is  always  less  than  one  for  typical  detonation 
combustion.  As  a  result,  the  typical  efficiency  of  a  detonation  cycle  is  greater  than  the 
efficiency  of  the  Brayton  cycle.  Based  on  initial  combustion  conditions  of  1  atmosphere 
and  300  °K,  the  Brayton  and  Humphrey  cycle  efficiencies  are  calculated  to  be  36.3  %  and 
56.4  %  respectively  with  a  difference  of  20.1  %  for  JP-lO/air  mixtures.  Thus 
thermodynamic  efficiency  of  a  PDE  combustion  cycle  is  very  high. 


D.  BASIC  ENGINE  OPERATION 

Cycle  frequency  is  defined  as  the  inverse  of  the  time  required  to  complete  a  full 
detonation  cycle.  Cycle  time  is  defined  as  the  time  required  to  fill  the  combustor  with 
fresh  fuel/air  mixture,  to  initiate  detonation,  to  allow  the  detonation  to  traverse  the  length 
of  the  tube  and  allow  the  unbumed  gases  to  expand  and  reduce  the  chamber  pressure  to 
permit  the  refill  process  to  occur  [Ref  3]. 
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(20) 


T  =T  +T  +T 

^  CYCLE  ^  DETONATION  ^  ^  EXPANSION  ^  ^  FILL 

The  instantaneous  peak  thrust  can  be  approximated  by  examining  the  pressure 
differential  acting  upon  the  interior  surfaces  of  the  combustor.  The  instantaneous  thrust 
[Ref  3]  which  acts  during  the  positive  P3,  or  head  end  pressure,  portion  of  the  cycle  can 
be  estimated  for  a  detonation  chamber  of  constant  cross-sectional  area,  A,  is  defined: 

FiNST=[_Pi{t)-Px'\^-Drag{t)  (21) 

Drag  forces  include  all  inviscid  and  viscous  forces  associated  with  flow  about  and 
through  the  vehicle.  The  average  thrust  [Ref  3]  is  defined:. 

(22) 

where 

Finst  =  instantaneous  or  unsteady  thrust 

f  =  cycle  frequency 

The  specific  impulse  (ISP)  can  be  computed  from  the  average  thrust  level  and  the 
fuel  consumption  or  mass  flow-rate  by  accounting  for  the  fuel  in  the  combustor  and 
second  by  accounting  for  the  fuel  and  auxiliary  oxygen  used  in  the  initiator  as  shown  in 
equation  (23)  [Ref  7]. 

V  =  7 - ^  (23) 

(wjy  +  -h  ni-o2_Aux  )S 

where 

f  =  cycle  frequency 

mf  =  mass  flow-rate  of  fuel  in  the  combustor 

nif  init  =  mass  flow-rate  of  fuel  in  the  initiator 

nio2_Aux  =  mass  flow-rate  of  auxiliary  oxygen  in  the  initiator 
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E.  DETONATION  INITIATION 

An  initiator  consists  of  a  small  auxiliary  combustor  filled  with  sensitive  mixture. 
The  unit  is  used  to  rapidly  generate  and  transmit  a  detonation  as  the  means  to  initiate  a 
detonation  in  a  larger  main  combustor  containing  less  sensitive  fuel/air  mixture.  Initiator 
concepts  vary  from  the  coaxial  to  transverse  or  splitter  a  plate  design  and  typically 
operate  on  fuel/oxygen  mixtures  or  a  blend  of  oxygen-enriched  air  as  the  oxidizer  [Ref 
7].  The  use  of  auxiliary  oxygen  provides  excellent  reliability,  repeatability  and  rapid 
ignition  but  its  minimization  is  required  since  it  is  treated  as  “fuel”  for  specific  impulse 
(Isp)  calculations  and  reduces  the  overall  system  performance  [Ref  7]. 

Some  examples  of  initiators  concepts  using  fuel/oxygen,  fuel/air  and  hybrid  are 
shown  in  Figure  13  [Ref.  8].  Concept  (a)  uses  an  oxygen/fuel  “plug”  as  the  head-end  of 
the  combustor;  which  works  well  but  uses  a  substantial  amount  of  oxygen.  Concepts  (b) 
and  (c)  use  a  fuel/oxygen  mixture  in  a  smaller  combustor  to  rapidly  generate  a  strong 
detonation  wave  to  initiate  the  fuel/air  mixture  in  the  combustor.  Concept  (b)  possesses  a 
solid  black  wall  at  the  diffraction  plane  whereas  concept  (c)  allows  the  wave  to  diffract 
initially  to  slight  larger  diameter.  The  splitter  plate  and  lateral  tube  concepts  shown  in  (d) 
and  (e)  respectively  have  been  explored  on  a  limited  scale.  Both  hybrid  designs  of 
concept  (f)  use  wall  shaping  and  shock  focusing  to  promote  the  generation  of  detonation 
wave.  Air  can  solely  be  used  as  the  oxidizer  in  both  designs  but  they  are  likely  to  be 
limited  in  the  operational  frequency,  reliability  and  their  requirement  for  additional  axial 
length  for  DDT  process  to  occur. 


Figure  13.  Examples  of  Initiator  Concepts. 
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Generally,  detonation  waves  produced  in  a  hydrocarbon  fuel/oxygen  initiator 
propagate  at  about  2800  to  3000  m/s  and  typically  near  1800  m/s  using  hydrocarbon 
fuel/air  mixtures.  Therefore  temperature  increases  across  the  detonation  shock  wave  for 
fuel/oxygen  initiators  are  on  the  order  of  3.5  times  more  than  for  fuel/air  initiators  based 
on  equations  (15)  and  (16). 

A  critical  area  of  concern  is  at  the  initiator  exit  plane  where  the  exiting 
detonation  wave  diffracts  into  the  combustor  and  is  the  motivation  for  characterizing  the 
effects  of  the  diffraction  condition  between  an  initiator  of  diameter,  Di,  and  the  main 
combustor  of  diameter,  D,  at  the  diffraction  plane.  Adequate  conditions  can  exist  through 
a  combination  of  diameter  ratio  (D/Di),  mixture  variation  (overdrive)  and  varying  degrees 
of  confinement  [Ref  7].  The  successful  and  unsuccessful  test  conditions  for  which 
detonation  transmission  occurred  for  the  ethylene-air  mixture  for  single-shot  cases  are 
shown  in  Figure  14  [Ref.  7].  It  shows  that  a  diffraction  ratio,  Y/Yq,  of  1.8  should  have  a 
successful  detonation  transition  from  the  initiator  to  the  combustor  for  an  initiator 
operating  at  equivalence  ratio  of  1.0. 


Figure  14.  Results  for  Two-Dimensional  Single-Shot  Diffraction  Tests. 
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The  critical  diameter  value  of  13  times  the  cell  size,  X,  of  a  mixture  has  been 
verified  many  times  for  the  transmission  of  a  detonation  wave  to  an  unconfmed  volume 
[Ref  7].  The  rule  holds  for  most  mixtures  and  is  especially  valid  for  mixtures  containing 
more  irregular  cell  spacing,  typically  fuel/air  mixtures  with  higher  activation  energies. 
Mixtures  containing  highly  regular  detonation  cell  structure,  such  as  argon  diluted 
fuel/oxygen  mixtures  often  require  a  larger  critical  diameter  than  the  13A,  rule,  thus 
emphasizing  the  importance  of  wave  front  structure  during  the  diffraction  processes  in 
generating  gas  dynamic  hot  spots  for  spontaneous  re-ignition  to  occur.  The  increased 
irregularity  in  the  cellular  structure  for  fuel/air  mixtures  often  assists  in  the  adjustment  to 
sudden  expansion  conditions  and  possessing  more  levels  of  instability,  therefore  more 
modes  by  which  spontaneous  re-initiation  may  occur  near  a  critical  diameter  value  [Ref. 

7]. 

Studies  on  the  re-initiation  mechanisms  of  shock  wave  reflections  at  a  rigid  wall 
from  the  propagation  of  a  quasi-detonation  in  an  obstacle-laden  channel  and  in  an 
imparting  spherical  blast  wall  on  a  rigid  wall  emphasize  the  importance  of  rapid  re¬ 
ignition  sites  immediately  behind  the  generated  Mach  stems  at  the  wall  [Ref.  7]. 
Additional  research  on  the  effects  of  overdriving  a  detonation  wave  during  diffraction 
from  a  smaller  combustion  tube  to  a  large  volume  have  shown  that  a  definitive  benefit 
exists  when  a  detonation  wave  is  allowed  to  propagate  into  a  less  reactive  mixture 
immediately  before  diffraction  occurs,  thus  creating  an  overdrive  condition  in  the  less 
reactive  mixture  [Ref  7].  Therefore  a  combination  of  shock  wave  reflections  and 
overdriven  conditions  would  be  the  mechanisms  that  govern  successful  initiator 
transmissions  on  most  PDEs. 

F.  VALVE-LESS  PULSE  DETONATION  ENGINE 

The  valve-less  PDE  developed  at  the  Naval  Postgraduate  School  uses  a 

continuous  air  flow  design,  which  eliminates  the  need  for  a  large  air  valve  to  supply  a 

fiiel/air  mixture  to  the  combustor  as  shown  in  Figure  15  [Ref  7].  The  initiator,  which  is 

similar  to  Figure  13(c),  operates  on  an  oxygen-enriched  fuel/air  mixture  to  rapidly  and 

reliably  generates  a  detonation  wave,  which  initiates  the  less  sensitive  fuel/air  mixture  in 

the  combustor.  Without  a  large  valve  on  the  air  flow,  a  convenient  flow  path  is  allowed  to 
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rapidly  fill,  detonate,  and  purge  the  combustion  chamber  at  rates  up  to  100  Hz  cycle 
frequency  [Ref.  7].  The  operational  cycle  of  the  valve-less  PDE  is  shown  below  in  Figure 
16  and  explained  in  Table  3  [Ref  7]. 


Figure  15.  Valve-less  PDE  Configuration. 


Figure  16.  Valve-less  PDE  Operational  Cycle 
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Stage 

Description 

A 

Air  flows  through  the  combustor  and  purge  the  previous  combustion  products 

B 

Fuel  is  injected  into  the  incoming  air  and  flows  into  the  combustor 

C 

Fuel/air  mixture  fully  fills  the  combustor 

D 

A  highly  detonable  mixture  is  rapidly  injected  into  the  initiator 

E 

Mixture  in  the  initiator  is  ignited  and  a  detonation  wave  forms 

F 

Detonation  wave  continue  to  propagate  in  the  initiator 

G 

Detonation  wave  exits  the  initiator  and  initiates  the  fuel/air  mixture  residing  in  the  combustor 

H 

After  the  detonation  exits  the  combustor,  a  series  of  rarefaction  waves  reduce  the  pressure  inside 
the  combustor  and  the  combustion  products  are  purged.  The  process  is  then  repeated. 

Table  3.  Valve-less  PDE  Operational  Cycle. 


Immediately  after  the  detonation  wave  from  the  initiator  diffracts  into  the 
combustor,  a  combustion-driven  shock  wave  begins  to  propagate  upstream  into  the 
incoming  air  stream.  The  propagation  of  the  shock  can  be  limited  to  a  certain  axial 
location  if  an  inlet  venturi  is  used.  After  the  detonation  wave  exits  the  main  combustor, 
rarefaction  waves  reduce  the  combustor  pressure  and  allow  the  purging  process  to  begin. 
The  process  then  repeats. 
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IV.  COMPUTATIONAL  FLUID  DYANMICS  TOOLS  AND 

TECHNIQUES 


A.  INTRODUCTION 

This  chapter  discusses  the  computational  fluid  dynamics  (CFD)  tools  used  and  the 
techniques  employed  in  modeling  the  various  combustors.  CFD  is  a  powerful  tool  for 
fluid  flow  analysis.  The  CFD  process  involves  many  primary  and  secondary  programs 
used  on  Silicon  Graphic  Iris  (SGI)  workstations  at  the  NFS.  The  interpretation  of  CFD 
results  requires  much  scientific  judgment  and  knowledge.  The  CFD  results  can  be 
improved  by  modifying  the  grid  mesh,  boundary  conditions,  numerical  differencing 
scheme,  time  step  and  more. 

B.  GRIDGEN  (GRID  GENERATION  PROGRAM) 

GRIDGEN,  version  9.6  is  a  grid  generation  program  developed  by  John  P. 
Steinbrenner  and  John  R.  Chawner  at  MDA  Engineering  and  supported  by  NASA  Ames 
Research  Center  via  Computer  Sciences  Corp  [Ref.  9].  It  is  a  software  system  for  the 
generation  of  3D,  multiple  block,  and  structure  grids  and  may  also  be  used  to  generate 
single  block  structured  grids,  single  surface  structure  grids  and  overset  structured  grids.  It 
is  an  interactive  code  used  to  decompose  a  three-dimensional  domain  into  blocks, 
distribute  grid  points  on  curves,  initialize  and  refine  grid  points  on  surfaces  and  initialize 
volume  grid  points.  It  is  menu  driven  but  does  not  have  advance  features  like  a  CAD 
system  to  model  complex  geometries.  However  it  can  read  IGES  files,  which  are  CAD 
generated. 


1.  Grid  Point  Distribution  Function 

A  hyperbolic  tangent  function  method  developed  by  Vinokur  is  based  on  a 
tangent  function,  which  either  has  a  read  (tangent)  or  a  purely  imaginary  (hyperbolic 
tangent)  argument.  The  method  allows  the  user  to  specify  both  the  spacing  between  the 
first  and  second  grid  points  and  the  spacing  between  the  next-to-last  and  last  grid  points. 
The  function  is  very  robust  and  results  in  smooth  transition  of  grid  point  clustering  along 

the  edge.  It  is  used  extensively  to  incrementally  distribute  the  grid  points  at  both  edges  of 
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the  connectors  along  the  axes  of  the  model  by  specifying  the  respective  distribution  factor 
at  each  edge  [Ref  9]. 

2.  Connectors 

Connectors  are  three-dimensional  space  curves,  which  represent  edges  of  the 
surface  grids  used  to  define  the  surfaces  or  domains.  A  connector  is  created  by  drawing 
the  shape,  assigning  the  number  of  grid  points  and  distributing  the  grid  points. 
Connectors  are  constructed  from  one  or  more  segments.  The  frequently  used  segments 
are  line,  three  point  circular  arc  and  conic  section  as  shown  in  Table  4  [Ref  9]. 


Type  of  Segment 

Description 

Line 

This  type  of  segment  passes  a  linear  polynomial  between  control  points.  At  least  two 
control  points  are  needed  to  specify  this  segment. 

Three  Point 
Circular  Arc 

This  type  of  segment  allows  the  user  to  define  a  circular  arc  of  not  more  than  1 80°  by 
selecting  three  control  points  that  lie  on  the  circle  as  shown  in  Figure  17  [Ref  9].  It  is 
constructed  that  begins  at  the  first  input  control  point,  ends  at  the  second  point,  and 
passes  through  a  third  control  point  elsewhere  in  space.  Cursor  movement  of  the  third 
control  point  is  constricted  so  that  circle  will  subtend  an  angle  of  not  more  than  180°. 

General  Conic 
Section 

This  type  of  segment  is  used  to  create  a  general  conic  section  spanning  an  arc  less 
than  180°  as  shown  in  Figure  18  [Ref  9].  It  is  constructed  from  three  (or  four)  input 
control  points.  The  resulting  conic  section  begins  and  ends  on  the  first  and  second 
control  point.  The  third  control  point  is  used  to  control  the  tangent  intersection  point, 
and  the  fourth  input  values  is  used  to  specific  the  p  value  that  determines  the  type  of 
conic  section. 

Table  4.  Type  of  Segment  and  Description. 


Figure  17.  Three  Point  Circular  Arc  Segment. 
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3.  Domains 

A  domain  represents  a  dimensionally  rectangular  patch  of  grid  points  which  is 
constructed  by  four  edges  with  each  edge  consists  of  a  string  of  connectors  that  connect 
end  to  end  as  shown  in  Figure  19  [Ref  9].  The  copy  menu  button  is  frequently  used  to 
copy  an  existing  domain  and  move  it  to  a  new  location. 


©  ©  © 


31 


4.  Blocks 

A  block  is  constructed  by  piecing  together  the  groups  of  domains  that  represent 
the  outer  shell  of  the  block.  It  consists  of  six  individual  faces  that  comprise  the  block 
shell  and  a  face  can  comprise  of  more  than  one  domain  as  shown  in  Figure  20  [Ref  9]. 


Figure  20.  Block  Construction  with  a  Face  with  Multiple  Domains. 

C.  GRIDED  (GRID  EDITING  PROGRAM) 

GRIDED,  version  3.0d  is  a  very  useful  menu-orientated  grid  editing  program.  It  is 
capable  of  converting  formatted  grid  files  to  unformatted  grid  files  and  vice  versa.  It  can 
covert  two  or  three-dimensional  whole  or  plane  grid  by  interchanging  the  grid  families, 
reserving  index  direction,  adding  extra  planes  by  reflected  symmetry  and  more.  An 
example  of  the  GRIDED  menu  is  shown  in  Appendix  D. 

D.  OVERFLOW  (FLOW  SOLVER  PROGRAM) 

The  OVERFLOW  program,  version  I.8w,  is  an  implicit  conservative  Navier- 
Strokes  (NS)  flow  solve.  It  is  a  complete  rewrite  of  F3D/Chimera  code  developed  by 
Joseph  E.  Steger  and  his  team  at  NASA  Ames  Research  Centre.  The  viscous  terms, 
which  are  partially  defined  by  the  shear  expressions,  are  set  to  zero  for  inviscid  test  cases. 
It  is  based  on  prefect  gas  law  algorithm,  thus  solution  will  be  stable  for  flow  less  than 
Mach  5.  It  is  not  capable  of  handling  chemical  reactions.  Table  5  shows  a  list  of  input  and 
output  files  commonly  used  and  generated  by  the  OVERFLOW  program  [Ref.  11]. 
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Filename 

File 

Type 

Description 

Filename.inp 

Input 

NAMELIST  inputs  and  boundary  conditions 

grid,  in 

Input 

Grid  file  in  single  or  multiple  grid  PLOT3D  format. 

q. restart 

Input 

Q  restart  file  in  modified  single  or  multiple  grid  PLOT3D  format. 

q.save 

Output 

Q  save  file  in  single  or  multiple  grid  modified  PLOT3D  format. 

q.bomb 

Output 

Q  bomb  file  in  modified  PLOT3D  format.  It  is  only  written  if  a  negative  density 
or  pressure  is  detected  or  if  the  specific  heat  ratio,  y,  falls  below  one. 

resid.out 

Output 

Flow  solver  residual  history. 

fomo.out 

Output 

Force,  moment  and  mass  flow  coefficient  history  per  wall. 

rpmin.out 

Output 

Minimum  density,  pressure  and  gamma  history. 

Filename.out 

Output 

Summary  of  NAMELIST  inputs,  boundary  conditions  and  status  of  solution 

Table  5.  OVERFLOW  Program  Input/Output  Files. 


The  input  parameter  specifications  to  OVERFLOW  program  are  done  using 
NAMELIST  inputs  as  shown  in  Appendix  A  and  are  written  in  the  input  files  such  as 
filename.inp.  Variables  generally  have  a  default  value  unless  specify  by  the  user.  The 
boundary  conditions  and  turbulence  model  specifications  for  the  input  files  are  shown  in 
Appendix  B.  A  batch  file  such  as  filename.batch  is  usually  used  to  run  the  OVERFLOW 
program  with  the  input  files  and  generate  the  output  files  such  as  filename. out  and  other 
output  files  listed  in  Table  5.  The  filename.out  output  file  contains  the  summary  of  the 
NAMELIST  inputs,  boundary  conditions  and  status  of  the  solution. 

The  various  flow  quantities  used  by  the  flow  solver  are  based  on  the  perfect  gas 
law  assumption  as  shown  below: 


p  =  pRT  =(y-l)pe.  =  {y-\)pcJ 

(24) 

=  yRT  =  ^ 

(25) 

p 

> 

II 

(26) 

c„ 

(27) 

R  =  Cp-c,={r-\)c^  =  ^ 

r 

(28) 
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Some  of  the  frequently  used  non-dimensional  thermodynamic  equations  [Ref.  11] 


are: 
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E.  FAST  (FLOW  SOLUTION  VISUALIZATION  PROGRAM) 

FAST  (Flow  Analysis  Software  Toolkit)  is  a  software  environment  for  data 
visualization.  A  series  of  module  runs  allows  the  user  to  examine  the  results  of  both 
numerical  and  experimental  simulation.  Users  can  load  data  files,  perform  data 
calculations,  visualize  the  calculated  results,  construct  scenes  of  three-dimensional 
objects,  and  plot,  animate  and  record  scenes.  A  list  of  frequently  used  FAST  modules  and 
their  functions  is  shown  in  Table  6. 


Modules 

Functions 

File  to 

It  loads  grid,  solution  and  PLOT3D  function  files  into  the  shared  memory  and  also  writes  out 
files.  Grids  can  be  structured  or  unstructured. 

Calculator 

It  is  used  to  calculate  CFD  functions.  PLOT3D  functions  are  offered,  and  custom  functions 
can  be  programmed.  Non-CFD  functions  can  also  be  created. 

Surfer 

It  visualize  the  grid  data  that  has  been  loaded  with  file  10  and  the  functions  that  have  been 
calculated  with  calculator  from  the  solution  data  loaded  with  file  10.  Surfer  renders  optional 
shading  using  interactive  system  light  sources. 

Viewer 

It  is  used  for  transformations  and  for  setting  lighting,  color,  clipping,  mirroring,  transparency, 
center  of  rotation,  anti-aliasing  and  scene  contents 

Plotter 

It  attaches  to  scalar  fields  in  shared  memory  or  reads  xy  ASCII  files  and  generates  two- 
dimensional  line  plots  and  histograms  in  computational  space. 

Table  6.  FAST  Modules. 
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F.  MODELING  TECHNIQUES 

The  approach  is  to  model  a  Mach  2.8  supersonic  shock  wave  transmitted  from  the 
initiator  to  the  combustor  through  an  initial  pressure  ratio  and  density  ratio  of  20  and  6 
respectively  between  the  initiator  and  combustor  using  the  OVERFLOW  flow  solver 
program.  A  total  of  ten  combustor  models  were  generated  with  four  cylindrical 
combustor  wall  geometries  and  six  scalloped  combustor  wall  geometries,  which  are 
described  later  in  the  chapter.  The  effects  of  shock  reflection/focusing  were  studied  by 
fixing  the  diameter  of  the  initiator  at  1.875  inches  and  varying  the  inner  diameter  of  the 
cylindrical  and  scalloped  combustor  at  2.25,  2.5,  2.75  and  3  inches,  forming  eight 
combustor  models.  Another  two  combustor  models  with  the  scallop  edges  of  the  2.75  and 
3  inches  scalloped  combustor  wall  geometries  modified  were  added. 


1.  Grid  Construction 

To  save  computational  space  and  time,  only  a  quarter  of  the  combustor  was 
modeled  using  the  GRIDGEN  program  with  both  the  extra  horizontal  and  vertical  planes 
being  defined  as  symmetry  planes  through  the  boundary  conditions  of  the  OVERFLOW 
program.  Figures  38  and  39  of  Appendix  C  show  the  grid  construction  models  of  a 
cylindrical  and  scallop  combustors  respectively.  Table  10  of  Appendix  C  shows  the 
coordinates  of  the  grid  points  of  the  combustors.  The  grid  mesh  constructions  of  both  the 
cylindrical  and  scallop  combustors  are  tabulated  in  Tables  11  and  12  of  Appendix  C. 
Both  the  general  conic  and  3  point  circular  arc  segments  were  used  to  construct  the 
quarter  arcs  and  the  scallop  arcs  of  the  combustors.  Figures  21  and  22  show  the  grid 
meshes  of  the  2.25  inches  diameter  cylindrical  and  scallop  combustors  respectively.  The 
red  and  blue  grid  mesh  represents  the  initiator  and  combustor  respectively. 

GRIDED  was  used  to  add  extra  planes  by  reflected  symmetry  on  the  horizontal 
and  vertical  plane  of  the  combustor,  increasing  the  grid  planes  in  the  circumferential 
index  (L  index)  from  81  to  83  so  that  the  OVERFLOW  program  could  recognize  the 
extra  planes  as  symmetry  planes  through  the  definition  of  the  boundary  conditions. 
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Figure  22.  Grid  Mesh  of  2.25  inches  Diameter  Scallop  Combustor. 


2.  Grid  Refinement 

Grid  refinement  is  more  like  an  art  than  a  science  with  no  definite  technique  or 
procedure.  The  Courant-Friedrichs-Levy  number,  CFL  is  defined  as: 
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(34) 
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Ax 


where 


At  =  time  step  of  the  integration,  known  as  DT  in  OVERFLOW  program 
Ax  =  grid  size 

a  =  wave  speed  of  the  following  first  order  wave  equation 

du  du  ^ 

—  +  a —  =  0 
dt  dx 


(35) 


The  convergence  and  stability  of  the  solution  depends  extensively  on  the 
combination  of  the  CFL  and  DT  and  finding  their  optimal  combination  may  differ  from 
case  to  case  depending  on  the  Mach  number.  The  general  guideline  is  that  grid  mesh  can 
be  refined  at  locations  of  sharp  edges,  sharp  comers,  flow  transition,  interest  and  both 
viscid  and  inviscid  walls. 


3.  Boundary  Conditions 

Table  13  of  Appendix  D  shows  the  boundary  conditions  of  the  combustors  used  in 
the  OVERFLOW  program.  The  core  axis  was  defined  as  axis  around  L  index  and  both 
the  extra  horizontal  and  vertical  grid  planes  were  defined  as  symmetry  planes.  The 
cylindrical  and  scallop  combustor  walls  were  defined  as  viscous  adiabatic  wall  with 
pressure  extrapolation  and  the  closed  rear  wall  of  the  initiator  was  defined  as  inviscid 
adiabatic  wall  with  pressure  extrapolation.  The  combustor  main  front  and  annular  rear 
outlets  were  defined  as  outflow  with  extrapolation.  The  boundary  conditions  used  were 
similar  to  the  validated  shock  tube  test  case  shown  in  the  README  file  of  Appendix  E. 


4.  Numerical  Differencing  Scheme 

Table  7  shows  a  list  of  numerical  differencing  schemes  used  in  the  validated 
shock  tube  test  case  as  shown  in  the  Appendix  G.  The  ability  of  Newton  sub-iterations  to 
decrease  factorization  errors  in  the  cross  flow  directions  was  tested.  Each  addition  of 
Newton  sub-iteration  took  the  same  amount  of  CPU  time  as  one  step  of  the  original 
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method.  Thus  a  ten  Newton  sub-iterations  run  would  take  ten  times  as  long  as  a  run  with 
no  sub-iteration. 


Base  Name 

No.  of  Newton 
Sub-Iterations 

Numerical  Differencing  Scheme 

Central 

Nil 

Central  Differencing,  Scalar  Dissipation,  Diagonal,  LHS 

Central  3 

3 

Central  Differencing,  Scalar  Dissipation,  Diagonal,  LHS 

Central  5 

5 

Central  Differencing,  Scalar  Dissipation,  Diagonal,  LHS 

Central  10 

10 

Central  Differencing,  Scalar  Dissipation,  Diagonal,  LHS 

Matrix 

Nil 

Central  Differencing,  Matrix  Dissipation,  Diagonal,  LHS 

Matrix  3 

3 

Central  Differencing,  Matrix  Dissipation,  Diagonal,  Left  Hand  Side 

Matrix  5 

5 

Central  Differencing,  Matrix  Dissipation,  Diagonal,  Left  Hand  Side 

Matrix  10 

10 

Central  Differencing,  Matrix  Dissipation,  Diagonal,  Left  Hand  Side 

Matrix  30 

30 

Central  Differencing,  Matrix  Dissipation,  Diagonal,  Left  Hand  Side 

Roe 

Nil 

Central  Differencing,  Matrix  Dissipation,  Diagonal,  Left  Hand  Side 

Roe  3 

3 

Roe  Upwinding,  Diagonal,  Left  Hand  Side 

Roe  5 

5 

Roe  Upwinding,  Diagonal,  Left  Hand  Side 

Roe  10 

10 

Roe  Upwinding,  Diagonal,  Left  Hand  Side 

Roe  30 

30 

Roe  Upwinding,  Diagonal,  Left  Hand  Side 

Lu  Roe  3 

3 

Roe  Upwinding,  LU-SGS  Algorithm,  Diagonal,  Left  Hand  Side 

Lu  Roe  10 

10 

Roe  Upwinding,  LU-SGS  Algorithm,  Diagonal,  Left  Hand  Side 

Table  7.  Numerical  Differencing  Scheme  of  Shock  Tube  Test  Case. 


The  density,  velocity  and  pressure  line  plots  of  the  various  scheme  of  the  test  case 
were  evaluated  against  the  experimental  results  of  Appendix  G.  It  was  observed  that  the 
central  differencing  scheme  coupled  with  scalar  dissipation  and  three  Newton  sub¬ 
iterations  was  sufficient  to  remove  factorization  errors.  Central  differencing  scheme  with 
matrix  dissipation  took  five  Newton  sub-iterations  and  Roe  upwinding  scheme  took  ten 
sub-iterations.  The  residual  files  showed  that  the  convergence  of  the  Newton  sub¬ 
iterations  slowed  down  greatly  after  the  first  three  iterations,  indicating  not  much  benefit 
from  additional  sub-iterations  for  the  larger  time  step  part  of  the  run.  However,  for  the 
first  ten  steps  with  the  smaller  DT,  additional  sub-iterations  continued  to  help 
convergence.  This  showed  that  the  time  step  of  0.01  might  be  too  large  for  maximum 
benefit  of  the  Newton  sub-iteration  procedure. 

The  numerical  differencing  scheme  chosen  for  all  the  combustor  models  was 
central  differencing  scheme  coupled  with  matrix  dissipation  and  ten  Newton  sub¬ 
iterations.  It  was  selected  based  on  its  test  case  solution  accuracy  against  the  experimental 
results  and  its  relatively  shorter  computation  time.  The  first  ten  steps  were  taken  with 
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time  step,  DT  of  0.0001  and  the  subsequent  steps  were  taken  with  time  step,  DT  of  0.001 
as  shown  in  central,  l.inp,  central.2.inp  and  central.S.inp  of  Appendix  F. 


5.  Input  Thermodynamic  Parameters  and  Assumptions 

All  the  combustor  models  had  the  same  initial  thermodynamic  parameters  with 
the  high  pressure,  temperature  and  density  region  filling  the  initial  half  of  the  initiator  and 
the  nominal  ambient  pressure,  temperature  and  density  region  filling  the  rest  of  the 
initiator  and  combustor  as  shown  in  Figure  23.  The  pressure  and  density  ratios  of  20  and 
6  between  these  regions  shown  in  Table  8  would  work  out  the  temperatures  to  be  960  °K 
and  288  °K  respectively  based  on  the  perfect  gas  relation.  The  specific  heat  ratios,  y  for 
both  regions  were  assumed  to  be  1.4  for  air.  These  initial  input  thermodynamic 
parameters  were  edited  in  the  makeq.f  file  as  shown  in  Appendix  F  and  was  later 
complied  to  makeq  executable  file  using  the  Fortran  77  complier.  The  execution  of  the 
makeq  file  would  generate  the  q.save  solution  file  with  the  initial  input  thermodynamic 
parameters  set  in  place.  The  central.  1  .batch  file  as  shown  in  Appendix  F  would  rename 
the  q.save  solution  file  to  q.restart  solution  input  file,  which  would  be  ready  to  feed  into 
the  OVERFLOW  flow  solver  program. 


High  Density,  Pressure  & 
Temperature  Region 

Nominal  Density,  Pressure 
&  Temperature  Region 

Density,  p  (kg/m^) 

7.350 

1.225 

Density,  p*  (non-dimension) 

6.00 

1.00 

Pressure,  p  (kPa) 

2026.500 

101.325 

Pressure,  p*(non-dimension) 

20.00 

1.00 

Temperature,  T  (°K) 

960.677 

288.150 

Specific  Heat  Ratio  of  Air,  y 

1.4 

1.4 

Table  8.  Initial  Input  Thermodynamic  Parameters 
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Figure  23.  Initial  Thermodynamic  Parameters  of  Combustors. 

6.  Flow  Solver  Program  Solutions  at  Different  Stages 

At  the  start  of  the  run,  the  high  density,  pressure  and  temperature  region  in  the 
initial  half  of  the  initiator  would  expand  and  generate  a  normal  shock  wave  propagating 
towards  the  opened  end  of  the  initiator.  The  pressure  behind  the  normal  shock  wave 
would  be  rapidly  reduced  to  a  plateau  pressure  of  approximately  six  atmospheres  due  to 
the  series  of  rarefaction  waves  that  forms  at  the  closed  end  of  the  initiator  as  shown  in 
Figure  24.  The  velocity  of  the  normal  shock  wave,  driven  by  the  reduced  plateau  pressure 
behind  it,  is  approximately  Mach  2.8  before  exiting  the  initiator. 

Subsequently,  the  normal  shock  wave  from  the  initiator  would  be  diffracted  into 
the  main  combustor.  The  diffracted  shock  would  be  reflected  by  the  combustor  wall 
geometry  and  would  form  two  normal  shock  waves,  one  propagating  upstream  into  the 
coming  air  stream  of  the  annular  combustor  outlet  and  another  propagating  downstream 
towards  the  main  outlet  of  the  combustor.  The  velocities  of  both  the  propagating  normal 
shock  waves  would  decrease  as  the  diameter  of  the  combustor  wall  increases. 

Due  to  the  limitations  of  the  flow  solver,  the  model  is  not  capable  of  simulating 
the  operating  condition  of  the  actual  PDF.  But  it  is  adequate  to  study  the  efforts  of 
reflected  shock/focusing  in  generating  regions  of  high  vorticity  and  subsequently  regions 
of  high  temperatures,  transmitted  from  the  initiator  into  the  combustor  through  a  weaker 
supersonic  shock  wave. 
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ConYerstion: 

P  =  P*xl.4atm  T=T*x143£23SK  p=p*xl.2kg^m3 


Length  (ft) 


Figure  24.  Thermodynamic  Parameters  Before  Normal  Shock  Wave  Exits  Initiator. 

7.  CFD  Procedure  Flow  Chart 

Figure  25  shows  a  CFD  procedures  flow  chart.  The  grid  mesh  was  generated 
using  GRIDGEN  and  edited  using  GRIDED.  Appendix  E  shows  the  procedures  of 
converting  grid.grd  to  grid.in  through  the  reserving  the  directions  of  the  grid  indexes, 
adding  extra  planes  by  reflected  symmetry  on  the  horizontal  and  vertical  plane  and  etc. 
FAST  was  used  to  check  the  grid  construction,  grid  distribution  and  directions  of  the  grid 
indexes.  The  executable  makeq  file  was  used  to  create  the  q.save  solution  file,  which  set 
up  the  initial  thermodynamic  parameters  and  the  q.save  solution  file  was  later  renamed  to 
q.restart  solution  input  file  through  the  execution  of  the  central.  1  .batch  file.  FAST  was 
used  to  check  the  initial  density,  pressure  and  temperature  values  of  the  combustor.  The 
high  density,  pressure  and  temperature  region  would  be  shaded  in  white  and  rest  in  black. 

There  were  three  input  files,  each  containing  identical  boundary  conditions.  The 
central. l.inp  input  file  had  10  steps  with  time  step,  DT,  of  0.0001  to  help  convergence. 
Both  the  central.2.inp  and  central.S.inp  input  files  had  the  same  time  step,  DT,  of  0.001 

but  a  different  number  of  steps,  NSTEPS,  of  ten  and  one  respectively  to  save  the 

41 


solutions  accordingly.  The  OVERFLOW  flow  solver  program  would  read  the  q.restart 
solution  input  file,  grid.in  grid  file  and  the  relevant  central. l.inp,  central.2.inp  or 
central.S.inp  input  file  and  initiate  the  run  through  the  execution  of  the  central.  1  .batch 
file.  The  solver  would  stop  the  run  at  each  designated  steps  as  specified  by  the 
central.  1  .batch  file,  save  and  copy  the  q.save  solution  file  to  a  new  file  name  to  prevent  it 
from  being  overwritten  and  rename  it  to  q.restart  solution  input  file.  The  program,  re¬ 
initiated  by  the  central.  1  .batch  file,  would  read  the  next  q.restart  solution  input  file, 
grid.in  grid  file  and  the  relevant  input  file  out  of  the  three  as  stated  above  and  the  cycle 
continues  until  the  last  command  of  the  central.  1  .batch  file  was  executed. 

If  the  solution  failed  due  to  a  negative  density  or  pressure  or  less  than  one  specific 
heat  ratio,  y,  the  previous  step  of  the  solution  would  be  saved  as  q.bomb  solution  file.  All 
the  central  output  files  and  the  residual  file  would  be  evaluated  to  determine  the  need  to 
increase  the  number  of  first  steps  and  reduce  the  time  step,  DT.  If  the  need  arises,  all  the 
central  input  files  would  be  modified  accordingly.  Subsequently,  FAST  was  used  to 
check  the  grid  mesh  locations  where  negative  density  or  pressure  or  the  less  than  one 
specific  heat  ratio,  y,  occurred.  That  region  of  the  grid  mesh  might  need  to  be  refined 
using  GRIDGEN.  The  process  of  determining  the  optimal  time  step,  DT,  and  grid  mesh 
size  was  indefinite  and  was  usually  based  on  trial  and  error  attempts. 
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Figure  25.  CDF  Procedures  Flow  Chart 
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V.  CFD  RESULTS  AND  ANALYSIS 


A.  INTRODUCTION 

This  chapter  discusses  and  analyzes  the  CFD  results  generated  by  the  flow  solver 
for  the  cylindrical  and  scalloped  combustor  geometries  described  earlier.  Interestingly,  it 
was  founded  that  the  maximum  temperature  regions  for  cylindrical  combustors  and 
smaller  diameter  scalloped  combustor  were  primarily  due  to  shock  focusing  efforts  as 
expected.  But  the  maximum  temperature  regions  of  the  larger  diameter  scalloped 
combustors  appeared  to  be  due  to  the  high  vorticity  regions  with  some  contributions  of 
shock  focusing  efforts. 


B.  TEMPERATURE 

1.  Duration  of  Maximum  Temperature  Regions 

An  annular  region  of  high  temperature  was  typically  formed  at  about  0. 1  to  0.2 
inches  after  the  initiator  and  between  the  outer  diameters  of  the  initiator  and  the 
combustor.  Figure  26  shows  the  maximum  temperature  of  one  of  the  high  temperature 
regions  and  its  duration  for  the  various  combustors. 

It  was  observed  that  for  the  cylindrical  combustors,  the  high  temperature  region 
would  increase  as  the  diameter  of  main  combustor  varied  from  2.25  to  2.75  inches,  but 
then  decreased  slightly  for  a  3  inches  diameter.  Generally,  the  maximum  temperature 
distributions  had  an  initial  sharp  rise,  followed  by  a  plateau  region  and  then  a  second 
gradual  rise  as  the  shock  wave  propagated  towards  the  end  of  the  main  combustor  except 
for  the  2.25  inches  diameter  scalloped  combustor.  The  2.25  inches  diameter  scalloped 
combustor  had  a  larger  second  rise  than  the  first.  The  high  temperature  regions  for  the 
second  temperature  rise  were  focused  at  the  center  of  each  scallop  of  the  combustor  wall 
geometry  as  shown  in  Figure  62  of  Appendix  J.  The  first  and  second  temperature  rises 
were  due  an  initial  weaker  shock  reflection,  followed  by  a  secondary  stronger  shock 
reflection. 

A  normal  shock  wave  would  propagate  at  Mach  2.8  in  the  initiator,  diffract  into 
the  combustor  and  propagate  at  Mach  2.5  there  after.  Based  on  the  initial  gas  temperature 
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and  normal  shock  properties,  the  temperature  across  the  normal  shock  propagating  at 
Mach  2.5  was  616°K.  Figure  26  shows  that  the  durations  for  the  maximum  temperature 
regions  defined  as  those  above  616°K,  were  approximately  50,  80,  160  and  40  micro¬ 
seconds  for  the  2.25,  2.5,  2.75  and  3  inch  diameter  cylindrical  combustors  respectively. 
The  temperature  duration  would  reach  its  optimal  at  2.75  inches  and  decline  drastically  at 
3  inches  and  thereafter  due  to  the  increasing  expansion  condition. 


Time  (us) 
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- D2.25m  Scalloped  - 02. Sin  Scalloped  - 02. 75in  Scalloped  - 03in  Scalloped 


Figure  26.  Maximum  Temperature  Vs  Duration  Plot  of  Various  Combustors 


The  scalloped  combustors  produce  maximum  temperature  durations  of 
approximately  120  and  40  micro-seconds  for  the  2.25  and  2.5  inches  diameter  geometries 
respectively.  The  maximum  temperature  duration  for  both  the  unmodified  and  modified 
scalloped  combustors  for  the  2.75  and  3  inches  diameters  could  not  be  accurately 
accounted  for  as  their  temperatures  rapidly  rise  beyond  4,000°K  preventing  their  CFD 
solutions  from  converging  due  to  the  inaccuracy  of  the  perfect  gas  law  at  those 
conditions. 
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The  Arrhenius  rate  law  of  equation  (13)  shows  the  relationship  of  reaction  rate 
with  temperature  rise.  The  maximum  temperature  regions  coupled  with  their  duration 
rapidly  accelerates  the  rate  of  chemical  reaction  of  a  combustible  mixture.  Therefore,  the 
2.75  and  3  inches  scalloped  combustors  would  have  considerable  potential  in  initiating 
the  less  sensitive  hydrocarbon  fuel/air  mixture  in  the  combustor  through  a  Mach  5.5 
shock  wave  from  the  hydrocarbon  fuel/air  initiator. 

Generally,  the  maximum  temperature  regions  of  the  scalloped  combustors  are 
approximately  1.5  to  5.5%  marginally  higher  than  the  cylindrical  combustor  for  the  2.25 
and  2.5  inches  diameter.  But  the  maximum  temperature  ratios  between  the  scalloped  and 
cylinder  combustors  are  2.4  and  more  than  3.2  for  the  2.75  and  3  inches  diameter 
respectively.  Thus  only  the  2.75  and  3  inches  diameter  cylindrical  and  scalloped 
combustors  are  selected  for  subsequent  comparisons. 


2.  Temperature  Variation  in  the  Main  Combustor 

The  temperature  plots  of  the  3  inches  diameter  cylindrical  and  scalloped 
combustors  are  compared  in  Figures  27  and  28.  The  temperature  variation  of  the  high 
temperature  regions  was  reasonably  uniform  for  the  cylinder  combustor  whereas  the 
extremely  high  temperature  regions  were  only  found  primarily  near  the  scallop  edges  of 
the  scalloped  combustors. 

Figure  29  shows  that  the  temperature  distributions  of  the  high  temperature  regions 
along  the  circumferential  index  (L  index)  of  the  grid  would  oscillate  circumferentially  for 
both  the  2.75  and  3  inches  diameter  cylindrical  combustors.  Whereas  for  the  2.75  and  3 
inches  diameter  scalloped  combustors,  extremely  high  temperature  spikes  (some  more 
than  4,000°K)  were  created  near  the  scallop  edges.  The  high  temperature  spikes  are  not 
exactly  consistent  at  each  scallop  edge  due  to  slight  dimensional  variation  of  each 
scallop.  The  majority  of  the  temperature  distributions  are  similar  in  magnitude  to  the 
cylindrical  combustor  results  for  both  scalloped  combustors.  Appendix  J  shows  the 
remaining  temperature  plots  of  the  other  combustors. 
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Figure  27.  Temperature  Plot  of  3  inehes  Diameter  Cylindrieal  Combustor. 
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Figure  28.  Temperature  Plot  of  3  inehes  Diameter  Sealloped  Combustor. 
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Figure  29.  Temperature  Variation  of  High  Temperature  Regions  for  2.75  and  3  inches 

Diameter  Cylindrical  and  Scalloped  Combustors. 


3.  Axial  Temperature  and  Pressure  Contour  Plots 

The  axial  temperature  and  pressure  contour  plots  of  the  3  inches  diameter 
cylindrical  and  scalloped  combustors  are  compared  in  Figures  30  and  31  at  203.15  and 
182.16  micro  seconds  after  the  initiation  of  the  run  respectively.  Appendix  1  showed  the 
remaining  temperature  and  pressure  contour  plots  of  the  other  combustors.  Although  not 
directly  shown,  after  the  shock  wave  from  the  initiator  expanded  into  the  combustor,  a 
portion  of  the  flow  diffracts  180°  and  forms  a  normal  shock  wave,  which  propagates 
towards  the  entrance  of  the  combustor. 

The  pressure  gradients  of  the  reflected  pressure  waves  were  similar  due  to  their 
constant  diameter  for  the  cylindrical  combustors,  but  the  reflected  shock  waves  from  the 
combustor  walls  were  not  ideal  since  the  incident  shock  waves  from  the  initiator  were 
highly  curves  and  not  sufficiently  planar  at  2.25  and  2.5  inches  diameter.  The  maximum 
reflected  shock  benefit  occurred  for  the  2.75  inches  diameter  combustor.  At  a  combustor 
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diameter  of  3  inches  and  larger,  this  magnitude  started  to  decrease  as  the  reflected 
pressure  region  increasingly  diffused  axially  to  the  rear  and  front  ends  of  the  combustor. 
The  pressure  contours  of  the  reflected  pressure  region  for  the  scalloped  combustors 
produced  interesting  results  along  the  scallop  edges;  which  represented  the  internal 
diameter  of  the  combustor  and  the  remainder  took  place  at  the  most  outer  tangents  of  the 
scalloped  wall. 

The  governing  mechanism  creating  the  high  temperature  regions  appeared  to  be 
primarily  shock  reflection  from  the  combustor  wall  and  not  vorticity  generation  for  the 
cylindrical  combustors  and  the  smaller  2.25  and  2.5  inches  diameter  scalloped 
combustors.  The  extremely  high  temperature  regions  of  the  2.75  and  3  inches  diameter 
scalloped  combustors  were  primarily  generated  by  high  vorticity  regions,  not  shock 
focusing  directly.  This  phenomenon  is  further  elaborated  later  in  this  chapter. 
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Figure  30.  Temperature  and  Pressure  Contour  Plot  of  3  inches  Diameter  Cylindrical 

Combustor. 
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Figure  3 1 .  Temperature  and  Pressure  Contour  Plot  of  3  inches  Diameter  Scalloped 

Combustor. 

C.  DURATION  OF  MAXIMUM  REFLECTED  PRESSURE  REGIONS 

Figure  32  shows  the  maximum  pressure  of  the  reflected  pressure  region  and  its 
duration  of  the  various  combustors.  Depending  on  the  diameter  and  geometry  of 
combustors,  the  maximum  reflected  pressure  distributions  would  have  an  initial  sharp  rise 
to  about  4.9  to  6  atmospheric  pressure  and  gradually  decline  to  2.5  to  3.4  atmospheres 
towards  the  end  of  the  shock  wave  propagation.  Generally,  the  scalloped  combustors  had 
slightly  higher  maximum  reflected  pressure  than  the  cylindrical  combustors  as  shown  in 
Figure  32. 
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Figure  32.  Maximum  Reflected  Pressure  Vs  Duration  Plot  of  Various  Combustors. 

Figure  33  shows  the  pressure  variation  at  the  high  temperature  regions  along  the 
circumferential  index  (L  index)  of  the  grid  for  the  2.75  and  3  inches  diameter  cylindrical 
and  scalloped  combustors.  Generally  the  pressures  of  the  scalloped  combustors  are  higher 
than  the  cylindrical  combustors,  but  not  significantly.  The  pressure  distributions  tend  to 
oscillate  circumferentially  as  well;  which  is  similar  to  their  temperature  distribution  but 
add  a  phase  shift  for  the  cylindrical  combustors.  The  erratic  pressure  distributions  were 
associated  to  their  corresponding  temperature  fluctuations  due  to  complexity  of  the 
scalloped  geometry. 


52 


Circ^Hinferential  \n4ex  (L  lndeK> 


-  -  -  -D3in  Cylindrical - D3in  Scalloped  -  -  -  -D2.75in  Cylindriical - D2.75in  Scalloped 


Figure  33.  Reflected  Pressure  Variation  Vs  Time  Plot  for  2.75  and  3  inches  Diameter 

Cylindrical  and  Scalloped  Combustors. 


D.  MACH  NUMBER 

Figures  34  and  35  show  the  Mach  number  vector  plots  of  the  3  inches  diameter 
cylindrical  and  scalloped  combustors.  Appendix  K  shows  the  Mach  number  plots  of  the 
other  combustors.  Both  cylindrical  and  scalloped  combustors  have  the  supersonic  Mach 
number  vectors,  expanded  outwardly  from  the  initiator,  and  the  sub-sonic  Mach  number 
vectors,  reflected  from  the  wall  of  the  combustor,  opposing  each  other. 

The  directions  of  both  opposing  Mach  number  vectors  shown  in  Figure  34  varied 
circumferentially,  revealing  large  counter-rotating  structures  for  the  3  inches  diameter 
cylindrical  combustor.  Figure  35  shows  that  the  Mach  number  vectors  crossed  multiple 
times  near  the  scallop  edges  for  the  3  inches  diameter  scalloped  combustor.  In  fact,  there 
are  three  opposing  Mach  number  vectors  near  the  scallop  edges  from  the  initiator,  the  left 
and  the  right  scallops  instead  of  two.  Thus  the  high  vorticity  generated  in  that  region  as 
described  later  in  this  chapter  was  responsible  for  causing  high  viscous  frictional  heat 
generation. 
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Figure  34.  Mach  No.  Vector  Plot  of  3  inches  Diameter  Cylindrical  Combustor 


Figure  35.  Mach  No.  Vector  Plot  of  3  inches  Diameter  Scalloped  Combustor 
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E.  VORTICITY  MAGNITUDE 

Figures  36  and  37  show  the  vorticity  magnitude  plots  of  the  3  inches  diameter 
cylindrical  and  scalloped  combustors.  The  vorticity  magnitude  of  the  cylindrical 
combustor  along  the  ring  of  the  sonic  region  is  uniformly  300  whereas  the  scalloped 
combustor  was  about  1100  near  the  scallop  edges.  Thus  the  vorticity  magnitude  ratio  of 
the  scalloped  over  the  cylindrical  combustors  was  about  3.67.  These  results  indicate  that 
although  shock  reflection  is  partially  responsible  for  producing  a  temperature  rise  at  the 
high  temperature  region,  the  vorticity  generation  and  subsequent  dissipation  mechanisms 
are  a  more  dominant  mechanism.  Appendix  L  shows  the  remaining  vorticity  magnitude 
plots  of  the  other  combustors. 
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Figure  36.  Vorticity  Magnitude  Plot  of  3  inches  Diameter  Cylindrical  Combustor. 
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Yorticity  Magnitude 


2000.0 

1800.0 

1600.0 

1400.0 

1200.0 

1000.0 

800.0 

600.0 

400.0 

200.0 

0.0 


Figure  37.  Vorticity  Magnitude  Plot  of  3  inches  Diameter  Scalloped  Combustor. 


F.  RESULTS  SUMMARY 

Table  9  shows  the  summary  of  results  for  the  10  combustors  evaluated.  The 
maximum  temperature  generated  by  the  scalloped  combustors  was  only  1.5  to  5.5% 
higher  than  the  cylindrical  combustor  for  the  2.25  and  2.5  inches  diameter  due  to  the 
effects  of  shock  reflection  and  reduction  in  vorticity  generation.  However,  the  maximum 
temperature  ratio  of  the  scalloped  over  the  cylindrical  combustors  were  2.4  and  more  than 
3.2  respectively  for  the  2.75  and  3  inches  diameter  due  to  the  large  amount  of  vorticity 
generated  and  the  supplemental  effects  of  shock  reflection/focusing.  Although  not 
validated,  prolonging  the  duration  of  these  high  temperature  regions  should  be  possible 
so  as  to  further  accelerate  the  chemical  reaction  of  the  combustible  mixture  in  those 
regions.  The  maximum  reflected  pressure  varied  about  minimally  from  a  nominal  value 
of  5.61  atmospheres  among  the  combustors. 
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_  Table  9.  Results  Summary  of  Combustors  _ _ 

Combustor  Time  Length  Temperature  Vorticity  Time  Length  Maximum 

(psec)  After  of  Hottest  Magnitude  (psec)  After  Pressure 
Initiator  Spot  (°K)  of  Hottest  Initiator  After 

(inch)  Spot  (inch)  Initiator 

_ (atm) 

2.25in  Diameter  159.84  O.I  745  350  150.49  0.9  5.74 

Cylindrical _ 

2.5in  Diameter  197.22  0.2  945  250  169.18  1.3  5.83 

Cylindrical _ 

2.75in  Diameter  205.64  0.2  1,336  300  178.53  1.4  5.45 

Cylindrical _ 

3in  Diameter  199.09  0.2  1,236  300  197.22  1.6  4.92 

Cylindrical _ 

2.25in  Diameter  205.64  0.1  885  300  150.49  0.9  5.77 

_ Scalloped _ 

2.5in  Diameter  171.99  0.2  917  600  150.49  0.9  5.96 

_ Scalloped _ 

2.75in  Diameter  163.57  0.2  3,202  1200  163.57  1.0  5.71 

Scalloped* _ 

3in  Diameter  178.79  0.2  >4,000  1100  178.53  1.1  6.01 

Scalloped* _ 

2.75in  Diameter  174.79  0.2  >4000  1750  169.18  1.2  5.88 

Modified  Scalloped* _ 

3in  Diameter  167.31  0.2  3,209  1100  167.31  1.0  4.84 

Modified  Scalloped* _ 

Notes: 

*  CFD  solution  fails  prematurely  due  to  extremely  high  temperature  generated,  which  defile  the  perfect  gas 
law  algorithm  of  the  OVERFLOW  flow  solver  program. 
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VI.  CONCLUSION  AND  RECOMMENDATION 


A.  CONCLUSION 

This  thesis  evaluated  the  effects  of  various  geometries  on  the  reflected  shock 
conditions  for  producing  regions  of  high  local  temperature  and  pressure.  A  total  of  10 
combustor  models  were  modeled  with  4  cylindrical  combustor  wall  geometries  and  6 
scalloped  combustor  wall  geometries.  The  effects  of  shock  refection/focusing  were 
studied  by  fixing  the  diameter  of  the  initiator  at  1.875  inches  and  varying  inner  diameter 
of  cylindrical  and  scalloped  combustor  at  2.25,  2.5,  2.75  and  3  inches,  forming  8 
combustor  models.  Another  2  combustor  models  with  the  edges  of  the  2.75  and  3  inches 
scalloped  combustor  wall  geometries  modified  were  added. 

The  results  revealed  that  the  reflected  pressure  and  high  temperature  regions  of 
both  the  cylindrical  and  scalloped  combustors  increase  as  the  combustor  diameter 
increases  eventually  to  a  point  of  diminishing  return.  The  maximum  reflected  pressure  of 
both  the  cylindrical  and  scalloped  combustors  did  not  vary  much  between  geometries. 
The  governing  mechanism  creating  the  high  temperature  regions  appears  to  be  primarily 
shock  reflection/focusing  from  the  combustor  wall  geometry  coupled  with  a  small 
amount  of  vorticity  generation  for  all  the  cylindrical  combustors.  This  phenomenon  is 
similar  for  the  smaller  2.25  and  2.5  inches  diameter  scalloped  combustors.  Whereas  the 
governing  mechanisms  of  creating  extremely  high  temperature  regions  for  the  2.75  and  3 
inches  diameter  scalloped  combustors  were  primarily  due  to  high  vorticity  generation 
near  the  scalloped  edges  which  supplemented  the  inherent  effects  of  shock 
reflection/focusing.  The  maximum  temperature  ratios  between  the  2.75  and  3  inches 
scalloped  and  cylinder  combustors  were  approximately  2.4  and  more  than  3.2 
respectively  due  to  a  large  vorticity  magnitude  ratio  of  more  than  3.67.  These  magnitudes 
should  be  viewed  as  approximate  values  since  the  CFD  flow  solver  program  utilizes  the 
ideal  gas  law  and  not  a  real  gas  model.  Therefore,  the  observed  threads  should  be 
consistent  but  reduced  in  magnitude  if  a  real  gas  model  is  utilized. 
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B.  RECOMMENDATIONS 

The  wall  contour  of  scalloped  combustor  should  be  modified  to  further  increase 
the  vorticity  magnitude,  particularly  the  angles  of  the  scallop  edges,  so  as  to  generate 
higher  temperature  and  duration.  A  flow  solver  program  capable  of  handling  chemical 
reaction  and  hypersonic  velocity  should  be  used  to  model  the  actual  PDE  conditions. 
Experimental  work  should  also  be  carried  out  to  verify  the  findings. 
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APPENDIX  A.  SUMMARY  OF  NAMELIST  INPUT  TO  OVERFLOW. 


SGLOBAL 

CHIMRA  TRUE-expect  grid  file  to  be  in  multiple  grid  format  and  PEGSUS-type  IBLANK 
and  interpolation  information  to  be  in  file  INTOUT  or  XINTOUT, 
FALSE-expect  single  grid  with  no  PEGSUS  information.  [FALSE] 

CDISC  TRUE-expect  to  read  a  NAMELIST  input  file  overdisc.coiiy  containing  CDISC 
inverse  design  control  information.  This  file  wUl  be  updated  by  OVERFLOW. 
FALSE-do  not  read  or  write  any  CDISC  information.  [FALSE] 

INCORE  TRUE-keep  data  for  all  grids  in  memory. 

FALSE-store  temporary  flowfield,  grid,  and  Chimera  interpolation  information  in 
scratch  files.  [FALSE] 

NSTEPS  Number  of  steps  to  advance  solution.  Use  zero  for  input  check.  [0] 

RESTRT  TRUE-read  restart  flowfield  from  file  q.restart. 

FALS E-start  from  initial  freestream  flowfield.  [FALSE] 

NSAVE  >0-Save  the  overall  solution  to  file  q.s&ve  every  how  many  steps. 

<  0-Save  the  solution  to  file  q.step#  every  how  many  steps.  [100] 

NFOMO  For  overset  force/moment  integration  package  (FOMOCO),  compute  integrated 
quantities  every  how  many  steps.  [10] 

NQT  Global  turbulence  model  type  declaration: 

0-algebraic  or  none. 

100-  1-equation  Baldwin-Barth  R'p  model. 

101- 1-equation  Spalart-Allmaras  R'p  model  (with  trip  line  specification). 

102- 1-equation  Spalart-Allmaras  R'p  model  (fully  turbulent). 

202- 2-equation  k-(jj  model. 

203- 2-equation  SST  model.  [0] 

NQC  Variable  y  model  type  declaration  (number  of  species): 

0-Constant  7,  1-gas  variable  7,  or  2-gas  variable  7  mixing  based  on  stagnation 
enthalpy. 

>2-Multiple  gas  variable  7  based  on  solution  of  NQC  species  continuity  equations. 

[0] 

MULTIG  Flag  to  enable/disable  multigrid  acceleration.  [FALSE] 

FMG  Flag  to  enable/disable  “full  multigrid”  (grid  sequencing)  precedure.  Only  for  use 
with  multigrid.  [FALSE] 

NGLVL  Number  of  multigrid  levels  to  use.  [3] 

FMGCYC(leveI#)  Number  of  steps  to  take  on  coarse  grid  levels  during  full  multigrid  procedure. 

(Here  index  1  is  the  coarsest  level,  2  the  next  finer,  etc.)  [300  for  all  coarse  levels] 
DTPHYS  Nondimensionalized  physical  timestep  for  time-accurate  simulations.  Nondimen- 
sionalization  is  based  on  grid  length  and  freestream  velocity.  [0.] 

NITNWT  Number  of  Newton/dual  subiterations  per  timestep,  or  0  for  no  subiteration.  [0] 
ORDNWT  (Not  used.  Number  of  orders  of  convergence  per  timestep  for  Newton/dual  subit¬ 
eration.)  [0.] 

FSONWT  1.0-first-order  in  time  Newton/dual  subiteration. 

2.0-secon<J-order  in  time  Newton/dual  subiteration. 

Intermediate  values  allowed.  [2.0] 
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$FLOINP 

ALPHA 

BETA 

FSMACH 

GAMLNF 

KEY 

PR 

PRT 

TINF 

RETINF 

XKINF 

FROUDE 

GVEC(3) 


TARGCL 

CLTARG 

CLALPH 


NTARG 


Angle-of-attack  (a),  deg.  [0.0] 

Sideslip  angle  (/3),  deg.  [0.0] 

Freestream  Mach  number  (Moo)*  [0.0] 

Freestream  ratio  of  specific  heats  (7  oo)*  [1*4] 

Reynolds  number  (Re)  (based  on  freestream  velocity  and  grid  length  unit).  [0.0] 
Prandtl  number  (Pt),  [0.72] 

Turbulent  Prandtl  number  (Prt).  [0.9] 

Freestream  temperature  (Too)j  deg.  Rankine.  [518.7] 

Freestream  turbulence  level  (Mt/Mi)oo  2-equation  turbulence  models.  [0.1] 

Freestream  turbulent  kinetic  energy  (^oo/V|o)  and  SST  turbulence  models. 

[0.0001] 

Ftoude  number  (gravity  term)  (Ri*)  (based  on  freestream  velocity  and  grid  length 
unit).  [0.0] 

Unit  up-vector  for  gravity  term.  (Note  that  this  vector  is  taken  verbatim-it  is  not 
modified  internaUy  by  the  angle-of-attack,  since  other  orientation  angles  (such  as 
bank  angle)  are  not  known.)  [0,0,1] 

TRUE-Enable  the  -driver  option.  [FALSE] 

Value  of  Ci^  which  the  code  will  try  to  match.  [0.] 

Initial  value  of  d(CL)/d( ALPHA)  used  to  estimate  the  necessary  change  in  AL¬ 
PHA.  This  is  only  used  for  the  first  correction,  after  which  the  CL  and  ALPHA 
changes  from  the  previous  correction  are  used.  [0.1] 

Number  of  steps  (or  cycles)  between  ALPHA  corrections,  with  the  following  ex¬ 
ceptions:  (1)  Corrections  are  not  done  during  grid  sequencing  (initial  coarse  grid 
cycles)  as  when  FMG=TRUE;  (2)  Corrections  are  not  done  on  the  first  or  last 
(fine  grid)  steps;  (3)  coarse  grid  cycles  are  not  counted  towards  NTARG.  [200] 
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$VARGAM 

IGAM  Options  for  specifying  calculation  of  7  when  not  solving  species  continuity  equa¬ 
tions  (i.e.,  NQC  <  2): 

0~Use  a  constant  7  value  of  GAMINF. 

1- Single  gas  with  variation  of  7  with  temperature  computed  using  ALTO-4,  AUTO- 
4. 

2- Two  gases,  variation  of  7  with  temperature  computed  using  ALTO-4,  AUTO-4; 
all  gas  1  below  HTl,  all  gas  2  above  HT2,  linear  mix  in  between.  [0] 

HTl  Total  enthalpy  ratio  which  the  mixture  is  all  gas  1.  [10.] 

HT2  Total  enthalpy  ratio  ^q/^Ooo  which  the  mixture  is  all  gas  2.  [10.] 

SCINF(gas#)  Freestream  species  mass  fraction  [1  for  gas  1,  0  for  all  others] 

SMVV(gas#)  Species  molecular  weight  MW,  or  normalized  molecular  weight  MWjMWoo  (if 
preferred).  [1.] 

ALT0(gas#)  Lower  temperature  range  polynomial  coefficient  ao  (540''R  <  T  <  ISOO^'R). 
[700/(700  -  1)] 

ALTl(gas#)  Lower  temperature  range  polynomial  coefficient  aj  (540® jR  <  T  <  1800®R).  [0.] 

ALT2(gas#)  Lower  temperature  range  polynomial  coefficient  (540®R  <  T  <  1800® jR).  [0.] 

ALT3(gas#)  Lower  temperature  range  polynomial  coefficient  (540® R  <T<  1800® H).  [0.] 

ALT4(gas#)  Lower  temperature  range  polynomial  coefficient  (540®  R  <  T  <  1800®/?).  [0.] 

AUT0(gas#)  Upper  temperature  range  polynomial  coefficient  uq  (1800®/?  <  T  <  9000®/?). 
[700/(700  -  1)] 

AUTl(gas#)  Upper  temperature  range  polynomial  coefficient  aj  (1800®/?  <  T  <  9000® /I).  [0.] 

AUT2(gas#)  Upper  temperature  range  polynomial  coefficient  02  (1800®/?  <  T  <  9000®/?).  [0.] 

AUT3(gas#)  Upper  temperature  range  polynomial  coefficient  03  (1800®/?  <  T  <  9000®/?).  [0.] 

AUT4(gas#)  Upper  temperature  range  polynomial  coefficient  04  (1800®/?  <  T  <  9000®/?).  [0.] 

SGRDNAM 

NAME  Grid  name  (for  input  identification  purposes  only,  not  used  within  the  code).  [“  ”] 
SNITERS 

ITER  >  0-Number  of  flow  solver  (sub)iterations  per  step  (number  of  iterations  before 
continuing  to  the  next  grid).  Each  iteration  of  the  flow  solver  implies  ITERT 
iterations  of  the  turbulence  model  and  ITERC  iterations  of  the  species  continuity 
equations. 

0-Do  not  advance  the  flow  solution;  take  ITERT  iterations  of  the  turbulence 
model,  ITERC  iterations  of  the  species  continuity  equations.  Use  ITER  =  ITERT  = 
ITERC  =  0  to  skip  this  grid.  [1] 
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SMETPRM 

IRHS  0-central  diifereiice  Euler  terms, 

2“Yee’s  symmetric  TVD  scheme, 

3- Liou  AUSM^-ap  upwind  scheme, 

4- Roe  upwind  scheme,  [0] 

ILHS  0-ARC3D  3'factor  block  tridiagonal  scheme. 

1- F3D  two“factor  scheme, 

2- ARC3D  3-factor  diagonal  scheme, 

3- LU-SGS  algorithm, 

4- D3ADI  algorithm  with  Huang  subiteration,  [2] 

ILHSIT  Number  of  subiterations  for  D3ADI  scheme,  [1] 

IDISS  1-F3D  dissipation  scheme  (2nd  and  4th  order  on  both  RHS  and  LHS  for  ILHS=2, 
2nd  order  only  on  LHS  for  ILHS=1). 

2- ARC3D  dissipation  scheme  (2nd  and  4th  order  on  both  RHS  and  LHS  for 
ILHS=2,  F3D  2nd  order  on  LHS  for  ILHS  =  1), 

3- TLNS3D  dissipation  scheme  (as  with  IDISS=2,  but  smooth  ph^  instead  of 

4- Matrix  dissipation  scheme  (see  additional  parameters  VEPSL,  VEPSN  and 
ROEAVG  in  SSMOACU),  [3] 

BIMIN  1, -Disable  low- Mach  number  preconditioning, 

<  1.— Enable  low-Mach  number  preconditioning. 

-*1, -Enable  low-Mach  number  preconditioning j  reset  BIMIN  to  3  x  [1.] 

MULTIG  Local  flag  to  en abie/disable  mnltigrid  acceleration,  [Default  is  value  of  MULTIG 
from  SGLOBAL.] 

SMOOP  Smoothing  coefficient  for  interpolation  (prolongation)  of  coarse  grid  solutions  onto 
the  next  finer  level  {FMG=TRUE  only),  [0,] 

SMOOC  Smoothing  coefficient  for  interpolation  of  multigrid  correction  onto  the  next  finer 
level,  [0,] 

SMOOR  Smoothing  coefficient  for  multigrid  residual  smoothing  before  restricting  to  the 
next  coarser  level,  [0,] 

CORSVI  Enable/disable  computation  of  viscous  terms  on  coarse  grid  levels.  [TRUE] 
$TIMACU 

ITIME  0-time  accurate, 

1 - local  timestep  scaling  (with  0,005  fudge  factor), 

2- local  timestep  scaling  (with  0  fudge  factor), 

3- Constanfc  CFL  number  (based  on  CFLMAX  value),  [1] 

ITIME=3  uses  the  sum  of  the  (max)  eigenvalue  in  each  coordinate  direction  to 
determine  the  local  CFL  number.  All  other  uses  of  CFLMIN/CFLMAX  use  the 
maximum  eigenvalue  to  determine  the  CFL  number, 

IRELAX  0— no  ovei  relaxation. 

1-N- relaxation  scheme  equable d  (Ref.  19),  [0] 

DT  Timestep  parameter.  Far  time-accurate  runs  (ITIME^O),  DT  is  nondimensional- 
ized  by  grid  length  and  freestream  speed- of- sound,  [0,5] 

TFOSO  1,0-first-order  time  accuracy  (Euler  implicit  scheme), 

2,0-second-order  time  accuracy  (trapeziodal  Implicit  scheme). 

Other  values  allowed;  0.5,  1.9  OK,  [1.0] 

CFLMIN  Minimum  CFL  number  cutoff,  [0.0] 

CFLMAX  Ma>dmum  CFL  number  cutoff.  [0,0  (none)] 
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SSMOACU 

ISPEC  Single  value  to  specify  ISPECJ,  ISPECK,  ISPECL.  Does  not  override  ISPECJ, 
ISPECK,  ISPECL.  [2] 

ISPECJ  — 1-sum  spectral  radii  in  J,K,L  for  dissipation  scaling  in  J. 

1- constant  coefficient  dissipation  in  J. 

2- spec tral  radius  in  J  for  dissipation  scaling  in  J. 

3- weighted  spectral  radius  in  J  (a  la  TLNS3D)  for  dissipation  scaling  in  J.  [ISPEC] 
ISPECK  — 1-sum  spectral  radii  in  J,K,L  for  dissipation  scaling  in  K. 

1- constant  coefficient  dissipation  in  K. 

2- spectral  radius  in  K  for  dissipation  scaling  in  K. 

3- weighted  spectral  radius  in  K  (a  la  TLNS3D)  for  dissipation  scaling  in  K.  flS- 
PEC] 

ISPECL  —  1-sum  spectral  radii  in  J,K,L  for  dissipation  scaling  in  L. 

1- constant  coefficient  dissipation  in  L. 

2- spectral  radius  in  L  for  dissipation  scaling  in  L. 

3- weighted  spectral  radius  in  L  (a  la  TLNS3D)  for  dissipation  scaling  in  L.  [ISPEC] 
SMOO  O.O-spectral  radius  is  computed  normally,  i.e.  as  (|I/|  -{-  /cc),  where 

^  =  y/ix  iy  Cly 

1.0-sound  speed  c  is  replaced  by  y  V^/Af|o  in  spectral  radius  terms. 
Intermediate  values  are  aUowed.  [1.0] 

DIS2  Single  value  to  specify  DIS2J,  DIS2K,  DIS2L.  Does  not  override  DIS2J,  DIS2K, 
DIS2L.  [2.0] 

DIS2J  2nd-order  smoothing  coefficient  in  J  (IDISS=2,3  only).  [DIS2] 

DIS2K  2nd-order  smoothing  coefficient  in  K  (IDISS=2,3  only).  [DIS2] 

DIS2L  2nd-order  smoothing  coefficient  in  L  (IDISS=2,3  only).  [DIS2] 

DIS4  Single  value  to  specify  DIS4J,  DIS4K,  DIS4L.  Does  not  override  DIS4J,  DIS4K, 
DIS4L.  [0.04] 

DIS4J  4th-order  smoothing  coefficient  in  J  (IDISS=2,3  only).  [DIS4] 

DIS4K  4th-order  smoothing  coefficient  in  K  (IDISS=2,3  only).  [DIS4] 

DIS4L  4th-order  smoothing  coefficient  in  L  (IDISS=2,3  only).  [DIS4] 

EPSE  Central  difference  smoothing  coefficient  (IDISS=1  only).  [0.05] 

EPSSGS  LU-SGS  left-hand  side  spectral  radius  epsilon  term  (ILHS=3  only).  [0.02] 

DELTA  Upwind  scheme  (MUSCL)  limiter  fix  parameter  (IRHS=3,4  only).  [1.] 

<  O.-Turn  off  limiter. 

O.-Koren  limiter. 

l.-Koren  limiter  with  CFL3D  parameter  €  =  S  (  • 

(Other  values  allowed.) 

FSO  1.0-first-order  for  flux  splitting/upwind  options. 

2 .0-second-order . 

3.0-third-order  (currently  implemented  only  for  Roe  and  AUSM"^-ap  upwinding). 
4.0-fourth-order  (implemented  only  for  central  differencing). 

Intermediate  values  allowed.  [2.0  for  IRHS=0,2,  3.0  for  IRHS=3,4] 

VEPSL  Matrix  dissipation  (minimum)  limit  on  linear  eigenvalues.  [0.] 

VEPSN  Matrix  dissipation  (minimum)  limit  on  nonlinear  eigenvalues.  [0.] 

ROEAVG  Matrix  dissipation  flag  for  using  Roe  averaging  for  half-grid  point  flow  quantities. 
[FALSE] 
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SVISINP 


Vise  TRUEi-include  aU  viscous  terms  including  cross  terms.  This  overrides  VISCJ, 
VISCK,  VISCL  and  VISCX. 

FALSE“include  only  specified  or  automatically  enabled  viscous  terms.  [FALSE] 
VISCJ  TRUE-include  viscous  thin  layer  terms  in  J. 

FALSE-viscous  terms  in  J  enabled  (only)  if  there  are  J-direction  viscous  walls. 
[FALSE] 

VISCK  TRUE-include  viscous  thin  layer  terms  in  K. 

FALSE-viscous  terms  in  K  enabled  (only)  if  there  are  Redirection  viscous  walls. 
[FALSE] 

VISCL  TRUE-include  viscous  thin  layer  terms  in  L. 

FALSE-viscous  terms  in  L  enabled  (only)  if  there  are  L-direction  viscous  walls. 
[FALSE] 

VISCX  TRUE-include  viscous  cross  terms  between  coordinate  directions  which  have  thin 
layer  terms  enabled. 

FALSE-no  viscous  cross  terms.  [FALSE] 

NTURB  Number  of  turbulent  wall  regions  to  be  specified  (optional).  [0] 

ITTYP(region#)  Type  of  turbulent  modeling  (see  Table  3.2). 

ITDIR(region#)  Coordinate  direction  of  turbulent  model  (away  from  wall  or  shear  layer).  1,  2,  3, 
—  1,  —2,  —3  represent  J,  K,  L,  — J,  —K,  — L,  respectively. 

JTLS(region#)  Starting  J  index. 

JTLE(region#)  Ending  J  index. 

KTLS(region#)  Starting  K  index. 

KTLE(region#)  Ending  K  index. 

LTLS(region^)  Starting  L  index. 

LTLE(region#)  Ending  L  index. 

TLPARl (region#)  Turbulence  model  parameter  (use  depends  on  turbulence  model  type). 

ITERT  Number  of  turbulence  model  (sub)iterations  per  flow  solver  iteration  (ITER),  or 
number  of  turbulence  model  (sub)iterations  if  ITER=0.  [1] 

ITLHIT  Number  of  subiterations  for  DDADI  scheme.  [1] 

CFLT  Turbulence  model  timestep  is  CFLT  times  the  flow  solver  timestep.  [l.j 

$BCINP 

NBC  Number  of  boundary  condition  regions  to  be  specified  (optional).  [0] 
IBTYP(region#)  Type  of  boundary  condition  (see  Table  3.3). 

IBDIR(region#)  Coordinate  direction  of  boundary  condition  (away  from  boundary  surface).  1,2, 
3,  —1,  —2,  —3  represent  J,  K,  L,  —J,  —K,  — L,  respectively. 

JBCS(region#)  Starting  J  index. 

JBCE(region#)  Ending  J  index. 

KBCS(region#)  Starting  K  index. 

KBCE(region#)  Ending  K  index. 

LBCS(region#)  Starting  L  index. 

LBCE(region#)  Ending  L  index. 

BCPARl(region#)  Boundary 'condition  parameter  (use  depends  on  boundary  condition  type). 
BCFILE(region#)  Filename  for  reading  boundary  data  (type  42). 

SSCEINP 

ITERC  Number  of  species  continuity  (sub)iterations  per  flow  solver  iteration  (ITER),  or 
number  of  species  continuity  (sub)iterations  if  ITER=0.  [1] 

CFLC  Species  convection  timestep  is  CFLC  times  the  flow  solver  timestep.  [1.] 
lUPC  0— 2nd-order  central  differencing  for  species  continuity  convection  terms. 

l-2nd-order  upwind  differencing  for  species  continuity  convection  terms,  [l] 
DIS2C  2nd-order  smoothing  coefficient  for  species  continuity  (JUPC  =  0  only).  [2.0] 
DIS4C  4th-order  smoothing  coefficient  for  species  continuity  {lUPC  =  0  only).  [0.04] 
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APPENDIX  B.  TURBULENCE  REGION  TYPES 


Typ  e  Des  cripti  on 

1  Bald  win- Lomax  boundary  layer  model  (with  variable  Degani-Schiif  cutoff) 
11  Baldwin-Lomax  shear  layer  model  (with  variable 

102  1’  or  2-equatioii  laminar  region  {zero  production) 

103  Spalart-Allmaras  boundary  layer  trip  line 


Type  Description 

1  In  viscid  adiabatic  wall  (pressure  extrapolation) 

2  Inviscid  adiabatic  wail  (normal  momentum  eqn) 

3  Inviscid  constant  temperature  wall  (pressure  extrapolation) 

4  Inviscid  constant  temperature  wall  (normal  momentum  eqn) 

5  Viscous  adiabatic  wall  (pressure  extrapolation) 

6  Viscous  adiabatic  wall  (normal  momentum  eqn) 

7  Viscous  constant  temperature  wall  (pressure  extrapolation) 

8  Viscous  constant  temperature  wall  (normal  momentum  eqn) 

10  Periodic  condition  (apply  to  either  1  or  last  plane) 

11  Symmetry  in  X  (apply  to  1  and/or  last  separately) 

12  Symmetry  in  Y  (apply  to  1  and/or  last  separately) 

13  Symmetry  in  Z  (apply  to  1  and/or  last  separately) 

14  Axis  (J  around) 

15  Axis  (K  around) 

16  Axis  {L  around) 

17  Symmetry  with  no  reflection  plane 

18  Periodic  flow /nonperiodic  grid  (apply  to  either  1  or  last  plane) 

21  2D  condition  in  Y  (3  planes  supplied,  ±1  in  Y) 

(apply  to  1  or  last  plane) 

22  Axisymmetric  condition  in  Y^  rotate  about  X  (3  planes  supplied, 

±1^^  rotation)  (apply  to  1  or  last  plane) 

30  Outflow  (extrapolation) 

31  Characteristic  condition  based  on  Riemann  invariants 

32  Supersonic/subsonic  inflow/outflow 

33  Specified  pressure  outflow 

35  Outflow  (Ist-order  extrapolation) 

40  Impose  freestream 

41  Nozzle  inflow  (poi  ^0  constant,  mass  flow  extrapolated) 

42  Prescribed  Q  (read  once  Irom  file) 

43  Prescribed  Q  with  siow-start  (read  from  file) 

44  Actuator  disk  (IDIR  is  flow  direction) 

45  Prescribed  Q  (read  from  file) /inflow- outflow  condition 

46  Prescribed  Q  with  slow  start  (read  from  file) /inflow-outflow  condition 

47  Freest  ream  /  c  harac  t  eris  tic  co  n  di  tion 

48  Simple  jet  mass  flow  condition 

49  Default  (no  change) 
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e  10.  Coordinates  of  Grid  Points  of  Combustors. 
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Table  1 1 .  Grid  Mesh  Construction  of  Cylindrical  Combustors. 


Face  Domain  Connector 

Type  of  Segment 

No.  of 
Grid  Points 

Distribution 

Function 

Distribution 
Factor  1,  ASl 

Distribution 
Factor  2,  ASl 

2.25in  /  2.50in  /  3.00in  Diameter  Cylindrical  Combustor 

A  A  a-a 

Line 

81 

Default 

Default 

Default 

a-b 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

a-c 

Line 

25 

Hyperbolic  Tangent 

O.lE-2 

O.lE-2 

b-c-d 

General  Conic  Section 

81 

Default 

Default 

Default 

B  b-c-d 

General  Conic  Section 

81 

Default 

Default 

Default 

b-e 

Line 

See  Notes 

Hyperbolic  Tangent 

O.lE-2 

O.lE-2 

c-1 

Line 

See  Notes 

Hyperbolic  Tangent 

0.1  E-2 

O.lE-2 

e-f-g 

General  Conic  Section 

81 

Default 

Default 

Default 

B  C  h-h 

Line 

81 

Default 

Default 

Default 

h-i 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

h-j 

Line 

25 

Hyperbolic  Tangent 

O.lE-2 

O.lE-2 

i-j-k 

General  Conic  Section 

81 

Default 

Default 

Default 

D  i-j-k 

General  Conic  Section 

81 

Default 

Default 

Default 

i-1 

Line 

See  Notes 

Hyperbolic  Tangent 

O.lE-2 

O.lE-2 

j-m 

Line 

See  Notes 

Hyperbolic  Tangent 

0.1  E-2 

O.lE-2 

1-m-n 

General  Conic  Section 

81 

Default 

Default 

Default 

C  E  a-a 

Line 

81 

Default 

Default 

Default 

h-h 

Line 

81 

Default 

Default 

Default 

a-h 

Line 

121 

Default 

Default 

Default 

a-h 

Line 

121 

Default 

Default 

Default 

D  F  a-b 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

h-i 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

a-h 

Line 

121 

Default 

Default 

Default 

b-i 

Line 

121 

Default 

Default 

Default 

G  b-e 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

i-1 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

b-i 

Line 

121 

Default 

Default 

Default 

e-1 

Line 

121 

Default 

Default 

Default 

E  H  a-c 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

h-j 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

a-h 

Line 

121 

Default 

Default 

Default 

c-j 

Line 

121 

Default 

Default 

Default 

I  c-f 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

j-m 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

c-j 

Line 

121 

Default 

Default 

Default 

f-m 

Line 

121 

Default 

Default 

Default 

F  J  e-f-g 

General  Conic  Section 

81 

Default 

Default 

Default 

1-m-n 

General  Conic  Section 

81 

Default 

Default 

Default 

e-1 

Line 

121 

Default 

Default 

Default 

f-m 

Line 

121 

Default 

Default 

Default 

Notes: 

No.  of  grid  points  for  2.25in  diameter  cylindrical  combustor 

is  11. 

No.  of  grid  points  for  2.50in  diameter  cylindrical  combustor  is  13. 

No.  of  grid  points  for  2.75in  diameter  cylindrical  combustor  is  16. 

No.  of  grid  points  for  3.00in  diameter  cylindrical  combustor  is  21 . 
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Table  12.  Grid  Mesh  Construction  of  Scallop  Combustors. 


Face  Domain  Connector 

Type  of  Segment 

No.  of 
Grid  Points 

Distribution 

Function 

Distribution 
Factor  1,  ASl 

Distribution 
Factor  2,  ASl 

2.25in  /  2.50in  /  3.00in  Diameter  Scallop  Combustor 

A  A  a-a 

Line 

81 

Default 

Default 

Default 

a-b 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

a-c 

Line 

25 

Hyperbolic  Tangent 

O.lE-2 

O.lE-2 

b-c-d 

General  Conic  Section 

81 

Default 

Default 

Default 

B  b-c-d 

General  Conic  Section 

81 

Default 

Default 

Default 

b-e 

Line 

See  Notes 

Hyperbolic  Tangent 

O.lE-2 

O.lE-2 

c-1 

Line 

See  Notes 

Hyperbolic  Tangent 

0.1  E-2 

O.lE-2 

e-f-g 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

f-h-i 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

h-j-k 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

j-l-m 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

B  C  n-n 

Line 

81 

Default 

Default 

Default 

n-o 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

n-p 

Line 

25 

Hyperbolic  Tangent 

O.lE-2 

O.lE-2 

o-p-q 

General  Conic  Section 

81 

Default 

Default 

Default 

D  o-p-q 

General  Conic  Section 

81 

Default 

Default 

Default 

o-r 

Line 

See  Notes 

Hyperbolic  Tangent 

O.lE-2 

O.lE-2 

p-y 

Line 

See  Notes 

Hyperbolic  Tangent 

O.lE-2 

0.1  E-2 

r-s-t 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

s-u-v 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

u-w-x 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

w-y-z 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

C  E  a-a 

Line 

81 

Default 

Default 

Default 

n-n 

Line 

81 

Default 

Default 

Default 

a-n 

Line 

121 

Default 

Default 

Default 

a-n 

Line 

121 

Default 

Default 

Default 

D  F  a-b 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

n-o 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

a-n 

Line 

121 

Default 

Default 

Default 

b-o 

Line 

121 

Default 

Default 

Default 

G  b-e 

Line 

See  Notes 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

o-r 

Line 

See  Notes 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

b-o 

Line 

121 

Default 

Default 

Default 

e-o 

Line 

121 

Default 

Default 

Default 

E  H  a-c 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

n-p 

Line 

25 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

a-n 

Line 

121 

Default 

Default 

Default 

c-p 

Line 

121 

Default 

Default 

Default 

I  c-1 

Line 

See  Notes 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

p-y 

Line 

See  Notes 

Hyperbolic  Tangent 

0.8E-2 

O.lE-2 

c-p 

Line 

121 

Default 

Default 

Default 

l-y 

Line 

121 

Default 

Default 

Default 

F  J  e-f-g 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

f-h-i 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

h-j-k 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

j-l-m 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

r-s-t 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

s-u-v 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

u-w-x 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

w-y-z 

3  Point  Circular  Arc 

21 

Default 

Default 

Default 

e-r 

Line 

121 

Default 

Default 

Default 

l-y 

Line 

121 

Default 

Default 

Default 

Notes: 

No.  of  grid  points  for  2.25in  diameter  scallop  combustor  is  1 

1. 

No.  of  grid  points  for  2.50in  diameter  scallop  combustor  is  13. 

No.  of  grid  points  for  2.75in  diameter  scallop  combustor  is  16. 

No.  of  grid  points  for  3.00in  diameter  scallop  combustor  is  2 

1. 
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APPENDIX  D.  BOUNDARY  CONDITIONS  OF  COMBUSTORS 


Table  13.  Boundary  Conditions  of  Cylindrical  and  Scallop  Combustors 


No. 

Location 

Boundary  Condition 
Type 

Boundary 

Condition 

Coordinate 

Direction 

Starting 

J 

Index 

Ending 

J 

Index 

Starting 

K 

Index 

Ending 

K 

Index 

Starting 

L 

Index 

Ending 

L 

Index 

1 

Core  Axis 

Axis,  L  Around 
(16) 

Positive  K 
(2) 

1 

(1) 

121 

(-1) 

1 

(1) 

1 

(1) 

I 

(1) 

83* 

(-1) 

2 

Initiator  Front 
Plane  Wall 

Inviscid  Adiabatic  Wall, 
Pressure  Extrapolation 
(1) 

Positive  J 
(1) 

1 

(1) 

1 

(1) 

1 

(1) 

25 

(25) 

1 

(1) 

83* 

(-1) 

3 

Combustor 
Front  Plane 

Outflow,  Extrapolation 
(30) 

Positive  J 
(1) 

1 

(1) 

1 

(1) 

26 

(26) 

See 

Notes 

(-1) 

1 

(1) 

83* 

(-1) 

4 

Combustor 

Horizontal 

Plane 

Symmetry  in  Z 
(13) 

Positive  L 

(3) 

1 

(1) 

121 

(-1) 

1 

(1) 

See 

Notes 

(-1) 

1 

(1) 

1 

(-1) 

5 

Combustor 
Vertical  Plane 

Symmetry  in  Y 
(12) 

Negative  L 

(-3) 

1 

(1) 

121 

(-1) 

1 

(1) 

See 

Notes 

(-1) 

I 

(1) 

1 

(1) 

6 

Combustor 
Rear  Plane 

Outflow,  Extrapolation 
(30) 

Positive  J 
(-1) 

121 

(-1) 

121 

(-1) 

1 

(1) 

See 

Notes 

(-1) 

1 

(1) 

83* 

(-1) 

7 

Combustor 

Inner 

Diameter  Wall 

Viscous  Adiabatic  Wall, 
Pressure  Extrapolation 
(5) 

Negative  K 
(-2) 

1 

(1) 

121 

(-1) 

See 

Notes 

(-1) 

See 

Notes 

(-1) 

1 

(1) 

83 

(-1) 

8 

Initiator  Inner 
Diameter  Wall 

Viscous  Adiabatic  Wall, 
Pressure  Extrapolation 
(5) 

Negative  K 
(-2) 

1 

(1) 

60 

(60) 

25 

(25) 

25 

(25) 

1 

(1) 

83 

(-1) 

9 

Initiator  Outer 
Diameter  Wall 

Viscous  Adiabatic  Wall, 
Pressure  Extrapolation 
(5) 

Positive  K 
(2) 

1 

(1) 

60 

(60) 

26 

(26) 

26 

(26) 

1 

(1) 

83 

(-1) 

Notes: 

The  ending  K  index  for  2.25in  diameter  cylindrical  and  scallop  combustors  is  35. 

The  ending  K  index  for  2.50in  diameter  cylindrical  and  scallop  combustors  is  37. 

The  ending  K  index  for  2.75in  diameter  cylindrical  and  scallop  combustors  is  45. 

The  ending  K  index  for  2.25in  diameter  cylindrical  and  scallop  combustors  is  50. 

Figure  40.  Boundary  Conditions  of  Combustors. 
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APPENDIX  E.  PROCEDURES  OF  CONVERTING  GRID.GRD  TO  GRID.IN. 

USING  GRIDED. 


******************************************************** 
GRIDED  version  3.0d 

******************************************************** 


Compile  time:  Fri  Mar  26  1 1:23:07  PST  1999 
Makefile  flag:  sgi 

Enter  input  PEOT3D  grid  filename:  grid..grd 

Enter  itin  (input  file  type)  and  itout  (output  file  type) 

1  =  3D  unformatted  whole,  2  =  3D  formatted  whole 
3  =  3D  unformatted  plane,  4  =  3D  formatted  plane 
5  =  2D  unformatted,  6  =  2D  formatted 
>2,1 

Iblanks  not  detected  in  grid  file. 

Current  grid  dimensions 

N  JMAX  KMAX  EMAX 
1  121  81  37 

Enter  operation  on  multiple  grid  file  from  following: 

0  =  Exit  from  this  menu  and  program 

1  =  Go  to  single  grid  menu 

2  =  Append  grids  by  reading  new  grid  file 

3  =  Remove  grids  from  current  list 

4  =  Substitute  a  grid  by  reading  new  single  grid 

5  =  Add/remove  iblanks  for  all  grids 

6  =  Swap  grid  numbers  of  two  grids 

7  =  Rearrange  sequence  of  all  grids 
9  =  Write  selected  grids  to  file 

1 

Enter  number  of  grids  to  perform  operations  on  (0  =  all) 
0 

Enter  iop  from  the  following  menu 
Operations  that  do  not  change  the  grid  dimensions 


1  =  Exit  from  this  menu  and  program 

2  =  Interchange  J  and  K  grid  families 

3  =  Interchange  K  and  E  grid  families 

4  =  Interchange  J  and  E  grid  families 

5  =  Reverse  index  direction  in  J  and/or  K  and/or  E 

6  =  Scale  followed  by  translate 

7  =  Mirror  about  x=0,  y=0  or  z=0  plane 

8  =  Rotate  about  x,  y,  z  or  any  arbitrary  axis 

9  =  Reset  axis/periodic/C-cut/constant-plane  bnd.  pts. 

Operations  that  change  the  grid  dimensions 


11=  Extract  a  subset 

12  =  Add  extra  planes  in  J,  K  or  E  by  extrapolation 

13  =  Add  extra  planes  in  J  or  K  by  reflected  symmetry 

14  =  Add  extra  planes  to  2D  grid  to  form  3D  grid  for 
'2D'  or  'axisymmetric'  option  in  OVERFEOW 

15  =  Duplicate  and  mirror  reflected  symmetry  grid  to  form  periodic  grid  in  K 

20  =  Concatenate  another  grid  to  current  grid 

21  =  Split  a  grid  into  2  grids  in  J,  K,  or  E  direction 

22  =  Create  a  new  grid  by  duplicating  a  grid  - 
further  ops  will  be  performed  on  duplicated  grid 

99  =  Go  to  multiple  grid  menu 
>5 

Enter  jrev,krev,lrev  (0=not  reverse,  l=reverse) 

>  1,0,0 

WARNING:  Options  2,3,4,5a,5b,5c,7  have  been  performed  an  odd  number  of  times. 
Grid  indices  have  changed  from  right-handed  to  left-handed  or  vice  versa. 

Do  you  wish  to  perform  further  operations  (l=yes,0=no)  >  1 
Enter  iop  from  the  following  menu 

Operations  that  do  not  change  the  grid  dimensions 


1  =  Exit  from  this  menu  and  program 

2  =  Interchange  J  and  K  grid  families 

3  =  Interchange  K  and  E  grid  families 

4  =  Interchange  J  and  E  grid  families 

5  =  Reverse  index  direction  in  J  and/or  K  and/or  E 

6  =  Scale  followed  by  translate 

7  =  Mirror  about  x=0,  y=0  or  z=0  plane 

8  =  Rotate  about  x,  y,  z  or  any  arbitrary  axis 

9  =  Reset  axis/periodic/C-cut/constant-plane  bnd.  pts. 


Operations  that  change  the  grid  dimensions 


11=  Extract  a  subset 

12  =  Add  extra  planes  in  J,  K  or  E  by  extrapolation 

13  =  Add  extra  planes  in  J  or  K  by  reflected  symmetry 

14  =  Add  extra  planes  to  2D  grid  to  form  3D  grid  for 
'2D'  or  'axisymmetric'  option  in  OVERFEOW 

15  =  Duplicate  and  mirror  reflected  symmetry  grid  to  form  periodic  grid  in  K 

20  =  Concatenate  another  grid  to  current  grid 

21  =  Split  a  grid  into  2  grids  in  J,  K,  or  E  direction 

22  =  Create  a  new  grid  by  duplicating  a  grid  - 
further  ops  will  be  performed  on  duplicated  grid 

99  =  Go  to  multiple  grid  menu 
>2 


Do  you  wish  to  perform  further  operations  (l=yes,0=no)  >  1 
Enter  iop  from  the  following  menu 

Operations  that  do  not  change  the  grid  dimensions 


1  =  Exit  from  this  menu  and  program 

2  =  Interchange  J  and  K  grid  families 

3  =  Interchange  K  and  E  grid  families 

4  =  Interchange  J  and  E  grid  families 

5  =  Reverse  index  direction  in  J  and/or  K  and/or  E 

6  =  Scale  followed  by  translate 

7  =  Mirror  about  x=0,  y=0  or  z=0  plane 

8  =  Rotate  about  x,  y,  z  or  any  arbitrary  axis 

9  =  Reset  axis/periodic/C-cut/constant-plane  bnd.  pts. 

Operations  that  change  the  grid  dimensions 


11=  Extract  a  subset 

12  =  Add  extra  planes  in  J,  K  or  E  by  extrapolation 

13  =  Add  extra  planes  in  J  or  K  by  reflected  symmetry 

14  =  Add  extra  planes  to  2D  grid  to  form  3D  grid  for 
'2D'  or  'axisymmetric'  option  in  OVERFEOW 

15  =  Duplicate  and  mirror  reflected  symmetry  grid  to  form  periodic  grid  in  K 

20  =  Concatenate  another  grid  to  current  grid 

21  =  Split  a  grid  into  2  grids  in  J,  K,  or  E  direction 

22  =  Create  a  new  grid  by  duplicating  a  grid  - 
further  ops  will  be  performed  on  duplicated  grid 

99  =  Go  to  multiple  grid  menu 

>  13 

Enter  jsyma,jsymb,ksyma,ksymb  (0=no  change,  l=add) 

>  1,0,0,0 

Enter  isymm  (l=x,2=y,3=z) 

>3 

Do  you  wish  to  perform  further  operations  (l=yes,0=no)  >  1 
Enter  iop  from  the  following  menu 

Operations  that  do  not  change  the  grid  dimensions 


1  =  Exit  from  this  menu  and  program 

2  =  Interchange  J  and  K  grid  families 

3  =  Interchange  K  and  E  grid  families 

4  =  Interchange  J  and  E  grid  families 

5  =  Reverse  index  direction  in  J  and/or  K  and/or  E 

6  =  Scale  followed  by  translate 

7  =  Mirror  about  x=0,  y=0  or  z=0  plane 

8  =  Rotate  about  x,  y,  z  or  any  arbitrary  axis 

9  =  Reset  axis/periodic/C-cut/constant-plane  bnd.  pts. 

Operations  that  change  the  grid  dimensions 


11=  Extract  a  subset 

12  =  Add  extra  planes  in  J,  K  or  E  by  extrapolation 

13  =  Add  extra  planes  in  J  or  K  by  reflected  symmetry 

14  =  Add  extra  planes  to  2D  grid  to  form  3D  grid  for 
'2D'  or  'axisymmetric'  option  in  OVERFEOW 

15  =  Duplicate  and  mirror  reflected  symmetry  grid  to  form  periodic  grid  in  K 

20  =  Concatenate  another  grid  to  current  grid 

21  =  Split  a  grid  into  2  grids  in  J,  K,  or  E  direction 

22  =  Create  a  new  grid  by  duplicating  a  grid  - 
further  ops  will  be  performed  on  duplicated  grid 

99  =  Go  to  multiple  grid  menu 
>5 
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Enter  jrev,krev,lrev  (O=not  reverse,  l=reverse) 

>  1,0,0 

WARNING:  Options  2,3,4,5a,5b,5c,7  have  been  performed  an  odd  number  of  times. 
Grid  indices  have  changed  from  right-handed  to  left-handed  or  vice  versa. 

Do  you  wish  to  perform  further  operations  (l=yes,0=no)  >  1 
Enter  iop  from  the  following  menu 

Operations  that  do  not  change  the  grid  dimensions 


1  =  Exit  from  this  menu  and  program 

2  =  Interchange  J  and  K  grid  families 

3  =  Interchange  K  and  E  grid  families 

4  =  Interchange  J  and  E  grid  families 

5  =  Reverse  index  direction  in  J  and/or  K  and/or  E 

6  =  Scale  followed  by  translate 

7  =  Mirror  about  x=0,  y=0  or  z=0  plane 

8  =  Rotate  about  x,  y,  z  or  any  arbitrary  axis 

9  =  Reset  axis/periodic/C-cut/constant-plane  bnd.  pts. 

Operations  that  change  the  grid  dimensions 


11=  Extract  a  subset 

12  =  Add  extra  planes  in  J,  K  or  E  by  extrapolation 

13  =  Add  extra  planes  in  J  or  K  by  reflected  symmetry 

14  =  Add  extra  planes  to  2D  grid  to  form  3D  grid  for 
'2D'  or  'axisymmetric'  option  in  OVEREEOW 

15  =  Duplicate  and  mirror  reflected  symmetry  grid  to  form  periodic  grid  in  K 

20  =  Concatenate  another  grid  to  current  grid 

21  =  Split  a  grid  into  2  grids  in  J,  K,  or  E  direction 

22  =  Create  a  new  grid  by  duplicating  a  grid  - 
further  ops  will  be  performed  on  duplicated  grid 

99  =  Go  to  multiple  grid  menu 

>  13 

Enter  jsyma,jsymb,ksyma,ksymb  (0=no  change,  l=add) 

>  1,0,0,0 

Enter  isymm  (l=x,2=y,3=z) 

>2 

WARNING:  Options  2,3,4,5a,5b,5c,7  have  been  performed  an  odd  number  of  times. 
Grid  indices  have  changed  from  right-handed  to  left-handed  or  vice  versa. 

Do  you  wish  to  perform  further  operations  (l=yes,0=no)  >  1 
Enter  iop  from  the  following  menu 

Operations  that  do  not  change  the  grid  dimensions 


1  =  Exit  from  this  menu  and  program 

2  =  Interchange  J  and  K  grid  families 

3  =  Interchange  K  and  E  grid  families 

4  =  Interchange  J  and  E  grid  families 

5  =  Reverse  index  direction  in  J  and/or  K  and/or  E 

6  =  Scale  followed  by  translate 

7  =  Mirror  about  x=0,  y=0  or  z=0  plane 

8  =  Rotate  about  x,  y,  z  or  any  arbitrary  axis 

9  =  Reset  axis/periodic/C-cut/constant-plane  bnd.  pts. 

Operations  that  change  the  grid  dimensions 


11=  Extract  a  subset 

12  =  Add  extra  planes  in  J,  K  or  E  by  extrapolation 

13  =  Add  extra  planes  in  J  or  K  by  reflected  symmetry 

14  =  Add  extra  planes  to  2D  grid  to  form  3D  grid  for 
'2D'  or  'axisymmetric'  option  in  OVEREEOW 

15  =  Duplicate  and  mirror  reflected  symmetry  grid  to  form  periodic  grid  in  K 

20  =  Concatenate  another  grid  to  current  grid 

21  =  Split  a  grid  into  2  grids  in  J,  K,  or  E  direction 

22  =  Create  a  new  grid  by  duplicating  a  grid  - 
further  ops  will  be  performed  on  duplicated  grid 

99  =  Go  to  multiple  grid  menu 
>5 

Enter  jrev,krev,lrev  (0=not  reverse,  l=reverse) 

>  1,0,0 _ 


Do  you  wish  to  perform  further  operations  (l=yes,0=no)  >  1 
Enter  iop  from  the  following  menu 

Operations  that  do  not  change  the  grid  dimensions 


1  =  Exit  from  this  menu  and  program 

2  =  Interchange  J  and  K  grid  families 

3  =  Interchange  K  and  E  grid  families 

4  =  Interchange  J  and  E  grid  families 

5  =  Reverse  index  direction  in  J  and/or  K  and/or  E 

6  =  Scale  followed  by  translate 

7  =  Mirror  about  x=0,  y=0  or  z=0  plane 

8  =  Rotate  about  x,  y,  z  or  any  arbitrary  axis 

9  =  Reset  axis/periodic/C-cut/constant-plane  bnd.  pts. 

Operations  that  change  the  grid  dimensions 


11=  Extract  a  subset 

12  =  Add  extra  planes  in  J,  K  or  E  by  extrapolation 

13  =  Add  extra  planes  in  J  or  K  by  reflected  symmetry 

14  =  Add  extra  planes  to  2D  grid  to  form  3D  grid  for 
'2D'  or  'axisymmetric'  option  in  OVEREEOW 

15  =  Duplicate  and  mirror  reflected  symmetry  grid  to  form  periodic  grid  in  K 

20  =  Concatenate  another  grid  to  current  grid 

21  =  Split  a  grid  into  2  grids  in  J,  K,  or  E  direction 

22  =  Create  a  new  grid  by  duplicating  a  grid  - 
further  ops  will  be  performed  on  duplicated  grid 

99  =  Go  to  multiple  grid  menu 
>2 

WARNING:  Options  2,3,4,5a,5b,5c,7  have  been  performed  an  odd  number  of  times. 
Grid  indices  have  changed  from  right-handed  to  left-handed  or  vice  versa. 

Do  you  wish  to  perform  further  operations  (l=yes,0=no)  >  1 
Enter  iop  from  the  following  menu 

Operations  that  do  not  change  the  grid  dimensions 


1  =  Exit  from  this  menu  and  program 

2  =  Interchange  J  and  K  grid  families 

3  =  Interchange  K  and  E  grid  families 

4  =  Interchange  J  and  E  grid  families 

5  =  Reverse  index  direction  in  J  and/or  K  and/or  E 

6  =  Scale  followed  by  translate 

7  =  Mirror  about  x=0,  y=0  or  z=0  plane 

8  =  Rotate  about  x,  y,  z  or  any  arbitrary  axis 

9  =  Reset  axis/periodic/C-cut/constant-plane  bnd.  pts. 

Operations  that  change  the  grid  dimensions 


11=  Extract  a  subset 

12  =  Add  extra  planes  in  J,  K  or  E  by  extrapolation 

13  =  Add  extra  planes  in  J  or  K  by  reflected  symmetry 

14  =  Add  extra  planes  to  2D  grid  to  form  3D  grid  for 
'2D'  or  'axisymmetric'  option  in  OVEREEOW 

15  =  Duplicate  and  mirror  reflected  symmetry  grid  to  form  periodic  grid  in  K 

20  =  Concatenate  another  grid  to  current  grid 

21  =  Split  a  grid  into  2  grids  in  J,  K,  or  E  direction 

22  =  Create  a  new  grid  by  duplicating  a  grid  - 
further  ops  will  be  performed  on  duplicated  grid 

99  =  Go  to  multiple  grid  menu 
>3 


Do  you  wish  to  perform  further  operations  (l=yes,0=no)  >  0 

Enter  output  PEOT3D  grid  filename:  grid.in 
Enter  0/1  for  single/multiple  grid  format 
0 

Output  grid  dimensions 

N  JMAX  KMAX  EMAX 
1  121  37  83 
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APPENDIX  F.  OVERFLOW  INPUT/OUTPUT  FILES 


CENTRAL.l. BATCH 

f77  -o2  -o  makeq  makeq.f 

makeq 

mv  q.save  q.restart 

/export/aahome/garth/overflowl.Sw.dir/overl.Sw/overflow  <  central. l.inp  >  central. Lout 
mv  q.save  q.restart 
cp  q.restart  q.save.  10a 

/export/aahome/garth/overflowl.Sw.dir/overl.Sw/overflow  <  central. 3.inp  >  central.3.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  10 

/export/aahome/garth/overflowl.Sw.dir/overl.Sw/overflow  <  central. 3.inp  >  central.3.out 
mv  q.save  q.restart 
cp  q.restart  q.save.20 

/export/aahome/garth/overflowl.Sw.dir/overl.Sw/overflow  <  central. 3.inp  >  central.3.out 
mv  q.save  q.restart 
cp  q.restart  q.save. 30 

/export/aahome/garth/overflowl.Sw.dir/overl.Sw/overflow  <  central. 3.inp  >  central.3.out 
mv  q.save  q.restart 
cp  q.restart  q.save.40 

/export/aahome/garth/overflowl.Sw.dir/overl.Sw/overflow  <  central. 3.inp  >  central.3.out 
mv  q.save  q.restart 
cp  q.restart  q.save. 50 

/export/aahome/garth/overflowl.Sw.dir/overl.Sw/overflow  <  central. 3.inp  >  central.3.out 
mv  q.save  q.restart 
cp  q.restart  q.save. 60 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central. 3.inp  >  central.3.out 
mv  q.save  q.restart 
cp  q.restart  q.save. 70 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central. 3.inp  >  central.3.out 
mv  q.save  q.restart 
cp  q.restart  q.save. 80 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central. 3.inp  >  central.3.out 
mv  q.save  q.restart 
cp  q.restart  q.save. 90 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central. 3.inp  >  central.3.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  100 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  101 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  102 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  103 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
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mv  q.save  q.restart 
cp  q.restart  q.save.  104 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  105 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  106 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  107 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  108 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  109 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  1 10 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  Ill 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  1 12 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 
cp  q.restart  q.save.  113 

/export/aahome/garth/overflowl.8w.dir/overl.8w/overflow  <  central.2.inp  >  central.2.out 
mv  q.save  q.restart 

cp  q.restart  q.save.XXX _ 


MAKEQ.F 


PARAMETER  ( jdim=121,kdim=42,ldim=83  ) 

o 

o 

II 

C 

endif 

C 

q(j,k,l,l)  =rinf 

real  q(jdim,kdim,ldim,6) 

q(j,k,l,2)  =  rinf^u 

real  dum(lO) 

q(j,k,l,3)  =rinf^v 

C 

q(j,k,l,4)  =  rinf^w 

C 

q(j,k,l,5)  =  pinf/gamma/(gamma-1.0)+0.5*(  q(j,k,l,2)**2+ 

gamma  =  1.4 

&  qG,k,l,3)**2+qG,k,l,4)**2  )/qG,k,l,l) 

C 

qG,k,l,6)  =  gamma 

C 

10  continue 

do  10j=l,jdim 

c 

do  10  k=l,kdim 

jmax  =  jdim 

do  10  l=l,ldim 

kmax  =  kdim 

if(j  .le.  30*jdim/121  .and.  k  .le.  25*kdim/42)  then 

Imax  =  Idim 

rinf  =  6.0 

open(unit=3,file='q.save',form='unformatted') 

pinf  =  20.0 

write(3)  jmax, kmax, Imax, 6,0 

u  =0.0 

write(3)  (dum(n),n=  1,10) 

V  =0.0 

write(3)  ((((qG,k,l,n),j= 1  ,jmax),k=  1  ,kmax),l=  1  ,lmax),n=  1 ,6) 

w  =  0.0 

close(3) 

else 

100  continue 

rinf=  1.0 

c 

pinf=  1.0 

stop 

u  =0.0 

end 

V  =0.0 
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CENTRAL.l.INP 


$GLOBAL 

SEND 

NSTEPS  =  10,  RESTRT  =  .T.,  NSAVE  =  10,  NITNWT  =  10, 

SSMOACU  SEND 

SEND 

SVISINP 

SFLOINP 

VISCK  =  .T.,  VISCJ  =  .T., 

FSMACH  =  0.5,  REY  =  9.E6,  TINF=  460.0, 

SEND 

SEND 

SBCINP 

SVARGAM  SEND 

NBC  =  9 

IBTYP  =  16,  1,  30,  13,  12,  30,  5,  5,  5, 

SGRDNAM 

IBDIR=  2,  1,  1,  3, -3, -1, -2, -2,  2, 

NAME  =  'Sod  shock  tube  grid'. 

JBCS  =  1,  1,  1,  1,  1,-1,  1,  1,  1, 

SEND 

JBCE  =-l,  1,  1,-1, -1,-1, -1,60,  60, 

SNITERS  SEND 

KBCS  =  1,  1,26,  1,  1,  1,-1,25,26, 

SMETPRM 

KBCE  =  1,25, -1,-1, -1,-1, -1,25,  26, 

IDISS  =  4, 

LBCS  =  1,  1,  1,  1,-1,  1,  1,  1,  1, 

SEND 

LBCE  =-l,-l,-l,  1,-1, -1,-1, -1,-1, 

STIMACU 

SEND 

ITIME  =  0,  DT  =  0.0001, 

SSCEINP  SEND 

CENTRAL.2.INP 


SGLOBAL 

SEND 

NSTEPS  =  1,  RESTRT  =  .T.,  NSAVE  =  1,  NITNWT  =  10, 

SSMOACU  SEND 

SEND 

SVISINP 

SFLOINP 

VISCK  =  .T.,  VISCJ  =  .T., 

FSMACH  =  0.5,  REY  =  9.E6,  TINF=  460.0, 

SEND 

SEND 

SBCINP 

SVARGAM  SEND 

NBC  =  9 

IBTYP  =  16,  1,  30,  13,  12,  30,  5,  5,  5, 

SGRDNAM 

IBDIR=  2,  1,  1,  3, -3, -1, -2, -2,  2, 

NAME  =  'Sod  shock  tube  grid'. 

JBCS  =  1,  1,  1,  1,  1,-1,  1,  1,  1, 

SEND 

JBCE  =-l,  1,  1,-1, -1,-1, -1,60,  60, 

SNITERS  SEND 

KBCS  =  1,  1,26,  1,  1,  1,-1,25,26, 

SMETPRM 

KBCE  =  1,25, -1,-1, -1,-1, -1,25,  26, 

IDISS  =  4, 

LBCS  =  1,  1,  1,  1,-1,  1,  1,  1,  1, 

SEND 

LBCE  =-l,-l,-l,  1,-1, -1,-1, -1,-1, 

STIMACU 

SEND 

ITIME  =  0,  DT  =  0.001, 

SSCEINP  SEND 

CENTRAL.3.INP 


SGLOBAL 

SEND 

NSTEPS  =  10,  RESTRT  =  .T.,  NSAVE  =  1,  NITNWT  =  10, 

SSMOACU  SEND 

SEND 

SVISINP 

SFLOINP 

VISCK  =  .T.,  VISCJ  =  .T., 

FSMACH  =  0.5,  REY  =  9.E6,  TINF=  460.0, 

SEND 

SEND 

SBCINP 

SVARGAM  SEND 

NBC  =  9 

IBTYP  =  16,  1,30,13,12,30,  5,  5,  5, 

SGRDNAM 

IBDIR=  2,  1,  1,  3, -3, -1, -2, -2,  2, 

NAME  =  'Sod  shock  tube  grid'. 

JBCS  =  1,  1,  1,  1,  1,-1,  1,  1,  1, 

SEND 

JBCE  =-l,  1,  1,-1, -1,-1, -1,60,  60, 

SNITERS  SEND 

KBCS  =  1,  1,26,  1,  1,  1,-1,25,26, 

SMETPRM 

KBCE  =  1,25, -1,-1, -1,-1, -1,25,  26, 

IDISS  =  4, 

LBCS  =  1,  1,  1,  1,-1,  1,  1,  1,  1, 

SEND 

LBCE  =-l,-l,-l,  1,-1, -1,-1, -1,-1, 

STIMACU 

SEND 

ITIME  =  0,  DT  =  0.001, 

SSCEINP  SEND 
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PEXT1.F 


C 

C  Extract  x  vs.  Cp  from  PLOT3D  grid  and  Q  files,  and  write  out  in 
C  XYPLOT  format. 

C  SET  FOR  MULTI  GRIDS 
C 

PARAMETER  (MI=20 1  ,MJ=80,MK=  1 23 ,MG=  I ) 

DIMENSION  XYZ(MI,MJ,MK, 3), Q(MI,MJ,MK, 5), ibl(mi,mj,nik) 

DIMENSION  JKL(MG,3) 

CHARACTER*80  FILEX,  FILEQ,  FILEO 
dimension  xp(5*MI),pp(5*MI) 

C 

DATA  GAMMA/I  .4/ 
data  xscale/I.O/ 

C 

pi  =  4.0*atan(  I.O  ) 

C 

C  Get  X  and  q  files. 

C 

WRITE(*,1) 

1  FORMATC  Enter  PLOT3D  grid  filename:  ',$) 

READ(*,2)  FILEX 

2  FORMAT(A) 

OPEN(UNIT=  I  ,FILE=FILEX,STATUS='OLD',FORM='UNFORMATTED') 
WRITE(*,4) 

4  FORMATC  Enter  PLOT3D  Q  filename:  ',$) 

READ(*,2)  FILEQ 

OPEN(UNIT=2,FILE=FILEQ,STATUS='OLD',FORM='UNFORMATTED') 

C 

C  Get  range  of  (i,j,k)  for  each  grid. 

C 

write(*,I2) 

12  formate  Enter  number  of  cuts  ',$) 
read(*,*)  ncuts 

do  1000  nc=I,ncuts 
WRITE(*,I4) 

14  FORMATC  Enter  XYPLOT  Cp  filename:  ',$) 

READ(*,2)  FILEO 

print  *, 'output  file  =  FILEO 

OPEN(UNIT=3,FILE=FILEO,STATUS='UNKNOWN',FORM='FORMATTED') 

write(*,I3) 

1 3  formate  Enter  number  of  grids  that  are  in  this  cut  ',$) 
read(*,*)  neg 

ii  =0 

do  500  nn  =  I,ncg 
WRITE(*,II) 

1 1  FORMATC  Enter  grid, IS, IE,  JS,JE,  KS,KE:  ',$) 

READ(*,*)  NG,IS,IE,JS,JE,KS,KE 
C 

C  Read  q  and  grid  file. 

C 

NGRIDS  =  1 

CRWT  READ(l)  NGRIDS 

READ(I)(JKL(N,I),JKL(N,2),JKL(N,3),N=I,NGRIDS) 

CRWT  READ(2)  NGRIDS 

READ(2)  (JKL(N,  I ), JKL(N,2), JKL(N,3),N=  I  ,NGRIDS) 
do  400  na=I,ngrids 
if(  na.eq.ng  )  then 
print  *,'na,ng  ',na,ng 
NI  =JKL(na,I) 

NJ  =JKL(na,2) 

NK  =JKL(na,3) 

IF  (NI.GT.MI  .OR.  NJ.GT.MJ  .OR.  NK.GT.MK)  THEN 
WRITE(*,3) 

3  FORMATCGrid  is  too  big,  exceeds  internal  dimensions.') 

GOTO  100 

ENDIF 

READ(1)((((XYZ(I,J,K,NX),I=1,NI),J=1,NJ),K=1,NK),NX=1,3) 

_ READ(2)  FSMACH,ALPHA,RE,TIME _ 
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CENTRAL.OUT 


OVERFLOW  -  OVERLAPPED  GRID  FLOW  SOLVER 
VERSION  1 .8w  12  December  2002 

Compiled  for  SINGLE  PRECISION:  using  _REAL 
Compile  time:  Thu  Feb  13  15:00:12  PST  2003 

Code  was  compiled  with  the  following: 

F77  =  f77 

F77FLAGS=  -mips4  -64  -03  -trapuv  -OPT:Olimit=0:IEEE_arithmetic=3:roundoff=3 
CC  =cc 

CFLAGS  =  -mips4  -64  -03  -trapuv 
CPP  =  /lib/cpp 
CPPFLAGS=  -DMIPS4 


Current  time:  Wed  Nov  19  21:23:31  2003 

GLOBAL  PARAMETERS  (SGLOBAL) 

CHIMERA  STYLE  INPUT?  (CHIMRA)  =  F 

RUN  INCORE?  (INCORE)  =  F 

RUNNING  CDISC  INVERSE  DESIGN?  (CDISC  )  =  F 

NUMBER  OF  STEPS  (NSTEPS)  =  1 0 

READ  RESTART  FILE?  (RESTRT)  =  T 

SAVE  RESTART  FILE  EVERY  (NSAVE  )  1 0  STEPS 

COMPUTE  FORCE/MOMENT  COEFS  EVERY  (NFOMO  )  1 0  STEPS 

TURBULENCE  MODEL  TYPE  (NQT)  =  0 

NUMBER  OF  SPECIES  (NQC)  =  0 

USE  MULTIGRID?  (MULTIG)  =  F 

USE  FULL  MULTIGRID?  (FMG)  =  F 

NO.  OF  GRID  LEVELS  (IF  MULTIG=.T.)  (NGLVL)  =  3 

NO.  OF  FMG  CYCLES  (IF  FMG=.T.)  (FMGCYC)  =00 

PHYSICAL  TIMESTEP  (NON-DIM  BY  V  INF)  (DTPHYS)  =  O.OOOOOE+00 

MAX  NUMBER  OF  NEWTON/DUAL  SUBITERS  (NITNWT)  =  1 0 

NO  ORD  CONVERGENCE  FOR  NEWTON/DUAL  SUB  (ORDNWT)  =  0.00000 

FIRST/SECOND  ORD  NEWTON/DUAL  SUB  (0-2)  (FSONWT)  =  2.00000 

ALLOCATING  MEMORY  FOR  GRID  AND  FLOWFIELD  ARRAYS: 


REQUESTING 

REQUESTING 

REQUESTING 

REQUESTING 


1 1388762  REAL  WORDS  FOR  FLOWFIELD  ARRAYS  (Q,S) 

6748896  REAL  WORDS  FOR  GRID  ARRAYS  (X,Y,Z, METRICS) 

421806  INTEGER  WORDS  FOR  GRID  ARRAY  (IBLANK) 

1 507660  REAL  WORDS  FOR  TEMPORARY  ARRAYS  (TMP,TMP2,TMP3) 


**  NOTE  **  Turning  off  force/moment  reporting  since  input  files 
(mixsur.fmp,  grid.ibi  and  grid.ptv)  do  not  exist. 

FLOW  CONDITIONS  ($FLOINP) 

ANGLE  OF  ATTACK  (DEG)  (ALPHA  )  =  0.00000 

SIDESLIP  (DEG)  (BETA  )  =  0.00000 

FREESTREAM  MACH  NUMBER  (FSMACH)  =  0.50000 

SPECIFIC  HEAT  RATIO  (GAMINE)  =  1 .40000 

REYNOLDS  NUMBER  (REY  )=  0.90000E+07 

PRANDTL  NUMBER  (PR  )=  0.72000 

TURBULENT  PRANDTL  NUMBER  (PRT  )=  0.90000 
FREESTREAM  TEMP  (DEG  R)  (TINF  )  =  460.00000 

VARIABLE  GAMMA  /  MULTIPLE  SPECIES  ($VARGAM) 
GAMMA  CALCULATION  METHOD  (0-2)  (IGAM  )  =  0 

TOTAL  ENTHALPY  RATIO  FOR  ALL  GAS  I  (HT I  )=  10.00000 
TOTAL  ENTHALPY  RATIO  FOR  ALL  GAS  2  (HT2  )=  10.00000 

INPUTS  FOR  GRID  I: 

GRID  NAME  (SGRDNAM) 

Sod  shock  tube  grid  (NAME  ) 

NUMBER  OF  TIME  STEPS  /  ITERATIONS  ($NITERS) 
ITERATIONS  PER  STEP  (ITER  )  =  I 

METHOD  CONTROL  PARAMETERS  ($METPRM) 

RIGHT-HAND-SIDE  OPTION  FLAG  (IRHS  )=  0 


82 


LEFT-HAND-SIDE  OPTION  FLAG  (ILHS  )=  2 

LEFT-HAND-SIDE  SUBITERATIONS  (ILHSIT)  =  I 
DISSIPATION  OPTION  FLAG  (IDISS  )  =  4 

LOW-MACH  PRECONDITIONING  PARAMETER  (BIMIN  )  =  1 .00000 

LOCAL  MULTIGRID  OPTION  (MULTIG)  =  F 

PROLONGATION  SMOOTHING  PARAMETERS  (SMOOPJ)  =  0.00000 

(IF  MULTIG=.T.)  (SMOOPK)  =  0.00000 

(SMOOPL)  =  0.00000 

CORRECTION  SMOOTHING  PARAMETERS  (SMOOCJ)  =  0.00000 

(IF  MULTIG=.T.)  (SMOOCK)  =  0.00000 

(SMOOCL)  =  0.00000 

RESIDUAL  SMOOTHING  PARAMETERS  (SMOORJ)  =  0.00000 

(IF  MULTIG=.T.)  (SMOORK)  =  0.00000 

(SMOORL)  =  0.00000 

USE  VISCOUS  TERMS  ON  COARSE  LEVELS?  (CORSVI)  =  T 
RECOMPUTE  MU  T  ON  FINEST  LEVEL?  (RECMUT)  =  F 

TIME  STEP/ACCURACY  PARAMETERS  ($TIMACU) 

TIME  STEP  SCALING  OPTION  FLAG  (0-2)  (ITIME  )  =  0 

RELAXATION  FACTOR  OPTION  FLAG  (O-I)  (IRELAX)  =  0 

TIME  STEP  (DT  )  =  O.IOOOOE-03 

FIRST/SECOND  ORDER  IN  TIME  (1-2)  (TFOSO  )  =  1.00000 

MINIMUM  CFL  NUMBER  (CFLMIN)  =  0.00000 

MAXIMUM  CFL  NUMBER  (CFLMAX)  =  0.00000 

SPATIAL  SMOOTHING/ACCURACY  PARAMETERS  ($SMOACU) 
SPECTRAL  RAD  SMOOTHING  (-1,1, 2, 3)  (ISPECJ)  =  2 

(ISPECK)  =  2 

(ISPECL)  =  2 

SPECTRAL  RAD  VS  VEL  SMOOTHING  (O-I)  (SMOO  )  =  1 .00000 

2ND  O  SMOOTHING  COEFS  (ARC3D-TYPE)  (DIS2J  )  =  2.00000 

(DIS2K)=  2.00000 
(DIS2L  )  =  2.00000 

4TH  O  SMOOTHING  COEFS  (ARC3D-TYPE)  (DIS4J )  =  0.04000 

(DIS4K )  =  0.04000 

(DIS4L  )  =  0.04000 

CENTRAL  DIFF  SMOOTHING  (F3D-TYPE)  (EPSE  )  =  0.05000 

LU-SGS  SPECTRAL  RADIUS  EPSILON  (EPSSGS)  =  0.02000 

ROE  LIMITER  FIX  PARAMETER  (DELTA  )  =  1 .00000 

CONVECTION  SPATIAL  DIFF  ORDER  (1-4)  (FSO  )  =  2.00000 

MATRIX  DISSIPATION  LINEAR  LIMIT  (VEPSL  )  =  0.00000 

MATRIX  DISSIPATION  NONLIN  LIMIT  (VEPSN )  =  0.00000 

USE  ROE  AVERAGE  IN  MATRIX  DISSIP?  (ROEAVG)  =  F 

VISCOUS/INVISCID  FLAGS  ($VISINP) 

INCLUDE  VISCOUS  TERMS  IN  J?  (VISCJ  )  =  T 
K?  (VISCK)  =  T 

L?  (VISCL  )  =  F 

INCLUDE  VISCOUS  CROSS  TERMS?  (VISCX  )  =  F 

TURBULENT  SURFACE  SPECS  ($VISINP) 

NUMBER  OF  TURBULENT  SURFACE  SPECS  (NTURB  )  =  -999 

TURBULENCE  MODEL  PARAMETERS  ($VISINP) 

ITERATIONS  PER  FLOW  SOLVER  ITERATION(ITERT  )  =  I 
LEFT-HAND-SIDE  SUBITERATIONS  (ITLHIT)  =  I 
LOCAL  TIME  STEP  CONSTANT  (CFLT  )  =  1 .00000 

BOUNDARY  CONDITION  SPECS  ($BCINP  ) 

NUMBER  OF  BOUNDARY  CONDITION  SPECS  (NBC  )=  9 

(IBTYP)  (IBDIR)  (JBCS)  (JBCE)  (KBCS)  (KBCE)  (LBCS)  (LBCE) 
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SPECIES  CONTINUITY  EQN  PARAMETERS  ($SCEINP) 

ITERATIONS  PER  FLOW  SOLVER  ITERATION(ITERC  )  =  I 
LOCAL  TIME  STEP  CONSTANT  (CFLC  )  =  1 .00000 

UPWIND  DIFFERENCING  FLAG  (O-I)  (lUPC  )  =  I 
2ND  O  SMOOTHING  COEF  (ARC3D-TYPE)  (DIS2C  )  =  2.00000 

4TH  O  SMOOTHING  COEF  (ARC3D-TYPE)  (DIS4C  )  =  0.04000 

GRID  SIZE  FOR  GRID  I : 

NUMBER  OF  POINTS  IN  J  (JD  )=  12 1 
K  (KD  )  =  42 
L  (ED  )=  83 

CHECKING  TIME  STEP  SPECIFICATION  FOR  GRID  I : 

RUNNING  TIME-ACCURATE  WITH  NEWTON  SUBITERATIONS 
WITH  DT  (BASED  ON  C  INF)  =  O.IOOOOE-03 

CHECKING  BOUNDARY  CONDITIONS  FOR  GRID  I : 

1)  BOUNDARY  CONDITION  TYPE#  16  DIRECTION  2 
AXIS  IN  K  (L  AROUND) 

1 .00  ORDER  EXTRAPOLATION  FOR  ZERO  SLOPE 
DIR=2  J-RANGE=  I  I2I  K-RANGE=  I  I  L-RANGE=  I  83 

2)  BOUNDARY  CONDITION  TYPE#  I  DIRECTION  I 
INVISCID  ADIABATIC  SOLID  WALL  (PRESSURE  EXTRAPOLATION) 

DIR=  I  J-RANGE=  I  I  K-RANGE=  I  25  L-RANGE=  I  83 

3)  BOUNDARY  CONDITION  TYPE#  30  DIRECTION  I 
OUTFLOW  (EXTRAPOLATION) 

DIR=  I  J-RANGE=  I  I  K-RANGE=  26  42  L-RANGE=  I  83 

4)  BOUNDARY  CONDITION  TYPE#  13  DIRECTION  3 
SYMMETRY  IN  Z,  IN  L 

DIR=3  J-RANGE=  I  12 1  K-RANGE=  I  42  L-RANGE=  I  I 

5)  BOUNDARY  CONDITION  TYPE#  12  DIRECTION  -3 
SYMMETRY  IN  Y,  IN  L 

DIR=-3  J-RANGE=  I  12 1  K-RANGE=  I  42  L-RANGE=  83  83 

6)  BOUNDARY  CONDITION  TYPE#  30  DIRECTION  -I 
OUTFLOW  (EXTRAPOLATION) 

DIR=-I  J-RANGE=  I2I  I2I  K-RANGE=  I  42  L-RANGE=  I  83 

7)  BOUNDARY  CONDITION  TYPE#  5  DIRECTION  -2 
VISCOUS  ADIABATIC  SOLID  WALL  (PRESSURE  EXTRAPOLATION) 

DIR=-2  J-RANGE=  I  I2I  K-RANGE=  42  42  L-RANGE=  I  83 

8)  BOUNDARY  CONDITION  TYPE#  5  DIRECTION  -2 
VISCOUS  ADIABATIC  SOLID  WALL  (PRESSURE  EXTRAPOLATION) 

NOTE  :  CONDITION  NOT  SPECIFIED  AT  PLANE  I  OR  LAST 
DIR=-2  J-RANGE=  I  60  K-RANGE=  25  25  L-RANGE=  I  83 

9)  BOUNDARY  CONDITION  TYPE#  5  DIRECTION  2 
VISCOUS  ADIABATIC  SOLID  WALL  (PRESSURE  EXTRAPOLATION) 

NOTE  :  CONDITION  NOT  SPECIFIED  AT  PLANE  I  OR  LAST 
DIR=2  J-RANGE=  I  60  K-RANGE=  26  26  L-RANGE=  I  83 

CHECKING  VISCOUS  AND  TURBULENCE  MODELING  SPECIFICATIONS  FOR  GRID  I : 

LAMINAR  OR  BALDWIN-LOMAX  O-EQUATION  TURBULENCE  MODEL  SELECTED 

INCLUDE  VISCOUS  TERMS  IN  J-DIRECTION 
INCLUDE  VISCOUS  TERMS  IN  K-DIRECTION 
NO  VISCOUS  TERMS  IN  L-DIRECTION 


FOR  GRID 

I 

AT  STEP 

10 

L2NORM  = 

0.7I552276E-05 

FOR  GRID 

I 

AT  STEP 

10 

L2NORM  = 

0.2I428687E-06 

FOR  GRID 

I 

AT  STEP 

10 

L2NORM  = 

0.50I77364E-08 

FOR  GRID 

I 

AT  STEP 

10 

L2NORM  = 

0.25I63888E-09 

FOR  GRID 

I 

AT  STEP 

10 

L2NORM  = 

0.455771 12E-I0 

FOR  GRID 

I 

AT  STEP 

10 

L2NORM  = 

0.32773270E-I0 

FOR  GRID 

I 

AT  STEP 

10 

L2NORM  = 

0.3065I485E-I0 

FOR  GRID 

I 

AT  STEP 

10 

L2NORM  = 

0.30269700E-I0 

FOR  GRID 

I 

AT  STEP 

10 

L2NORM  = 

0.30I996I4E-I0 

FOR  GRID 

I 

AT  STEP 

10 

L2NORM  = 

0.30I89667E-I0 

FLOW  SOLVER/TURB  MODEL/SPECIES  CONVERGENCE  RATES :  1 .0889  0.0000  0.0000 

Current  time:  Wed  Nov  19  21:39:16  2003 
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APPENDIX  G.  README  FILE  OF  TEST  CASE. 


This  is  a  shock  tube  test  case  for  time  accuracy. 

Problem  setup  (grid,  initial  conditions  and  exact  solution)  comes 
courtesy  of  Mr.  Robert  Tramel,  Micro  Craft,  Inc.,  Tullahoma,  TN. 

The  problem  is  set  up  as  a  3-dimensional  problem  (grid  is  101x51x5), 
with  inviscid  sidewalls  and  viscous  walls  top  and  bottom.  This  tests 
the  ability  of  the  Newton  subiteration  to  drive  out  factorization  errors 
in  the  crossflow  directions.  Note  that  each  additional  Newton 
subiteration  (beyond  1)  takes  the  same  amount  of  CPU  time  as  one  step  of 
the  (original)  method.  Thus  running  for  10  sub  iterations  takes  10  times 
as  long  as  running  without  subiterations. 

Ten  iterations  are  taken  with  a  timestep  (DT)  of  0.001  to  start  the 
simulation,  then  18  steps  are  taken  with  a  timestep  of  0.01,  for  a  total 
(nondimensional)  time  advance  of  0. 19.  (Time  is  nondimensinalized  using 
the  freestream  speed  of  sound  and  one  unit  of  length  in  the  CFD  grid.) 

The  initial  smaller  timesteps 

are  needed  (at  least  for  central/scalar)  to  keep  the  solution  from 
bombing  due  to  the  initial  transient. 

All  runs  are  done  with  2nd-order  (3 -point  backward)  time  accuracy, 
except  for  initial  steps  of  each  run  (and  when  DT  is  changed).  The 
previous  time  level  (Q^(n-l))  is  NOT  stored  in  the  q.save/q.restart  file. 


Running  cases: 

Cases  can  be  run  using  the  "runcase"  script  (and  results  for  a  particular 
case  deleted  with  "cleancase").  Scripts  take  a  single  argument  which  is 
the  base  name  for  two  input  files,  (basename).l.inp  and 
(basename).2.inp.  Results  are  (1)  an  XYPLOT  file  of  x  vs.  rho,  V  and  p 
in  shock. (basename),  (2)  a  residual  history  in  resid.(basename),  (3)  a  Q 
file  in  q.(basename),  and  (4)  a  PLOT3D  command  file  in  (basename). com 
for  plotting  density  contours. 

The  script  "runall"  runs  a  selection  of  the  following  cases;  "plotall" 
generates  some  comparison  plots  (using  XYPLOT  and  PLOT3D);  and 
"cleanall"  cleans  up  these  cases. 

Input  for  the  following  cases  is  included: 


basename  subiter 

differencing  method 

central 

no 

central,  scalar  diss,  diagonal  LHS 

centraB 

3 

central,  scalar  diss,  diagonal  LHS 

centrab 

5 

central,  scalar  diss,  diagonal  LHS 

central  10 

10 

central,  scalar  diss,  diagonal  LHS 

matrix 

no 

central,  matrix  diss,  diagonal  LHS 

matrix3 

3 

central,  matrix  diss,  diagonal  LHS 

matrix5 

5 

central,  matrix  diss,  diagonal  LHS 

matrix  10 

10 

central,  matrix  diss,  diagonal  LHS 

matrix30 

30 

central,  matrix  diss,  diagonal  LHS 

roe 

no  Roe,  diagonal  LHS 

roe3 

3  Roe,  diagonal  LHS 

roe5 

5  Roe,  diagonal  LHS 

roe  10 

10  Roe,  diagonal  LHS 

roe30 

30  Roe,  diagonal  LHS 

lu_roe3 

3 

Roe,  LU-SGS  LHS 

luroelO 

10 

Roe,  LU-SGS  LHS 

Plotting  results: 

The  following  PostScript  files  have  already  been  prepared: 
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method.ps  Compares  central  differencing,  central  with  matrix 

dissipation,  and  Roe  up  winding  (cases  centraB,  matrix5, 
roe  10). 

central.ps  Compares  different  numbers  of  subiterations 
(cases  central,  centraB,  centraB). 
matrix.ps  Compares  different  numbers  of  subiterations 
(cases  matrix3,  matrix5,  matrix  10). 
roe.ps  Compares  different  numbers  of  subiterations 
(cases  roe5,  roe  10,  roe30). 

roe_lhs.ps  Compares  Roe  upwind  with  different  left-hand  sides 

(cases  roe  10,  lu  roelO). 

If  you  have  XYPLOT,  results  can  be  compared  by  typing 

xyplot  -f  shock.fmt  shock.exact  shock.(casel)  shock.(case2) ... 

where  (easel),  (case2)  are  to  be  replaced  by  the  desired  basenames. 

Density  contours  in  the  3D  grid  can  be  plotted  using  PLOT3D,  with  the 
(basename).com  command  file.  This  will  show  the  deviation  of  the 
solution  in  the  crossflow  directions. 


Some  observations: 

o  For  standard  central  differencing  (scalar  4/2  dissipation), 

3  subiterations  is  adequate  to  remove  factorization  errors  (as 
determined  by  looking  at  the  line  plots  of  density,  velocity  and 
pressure,  and  by  looking  at  contour  plots  of  density  in  K=26  and  L=3 
planes). 

o  Similarly,  central  differencing  plus  matrix  dissipation  takes  5 
subiterations  (with  slight  changes  past  this),  and  Roe  upwinding  takes 
10  subiterations. 

o  At  10  sub  iterations,  the  (diagonalized)  central  left-hand  side 
(ILHS=2)  appears  to  have  less  crossflow  distortion  than  the  LU-SGS 
scheme  (ILHS=3).  Specifically,  the  contours  near  the  viscous  walls 
are  much  straighter. 

o  Looking  at  the  residual  files,  convergence  of  the  Newton  subiterations 
slows  down  greatly  after  the  first  3  iterations,  indicating  not  much 
benefit  from  additional  subiterations  (for  the  larger  timestep  part  of 
the  run).  However,  for  the  first  10  steps  with  the  smaller  DT, 
additional  subiterations  continues  to  help  convergence.  This 
indicates  that  the  timestep  of  0.01  may  simply  be  too  large  for 
maximum  benefit  of  the  Newton  subiteration  procedure.  An  examination 
of  accuracy  vs.  timestep  is  in  order. 


Reference  for  the  physical  problem: 

G.A.  Sod,  "A  Survey  of  Several  Finite  Difference  Methods  for  Systems 
of  Non-Linear  Hyperbolic  Conservation  Laws,"  JCP,  Vol.  27,  1978, 
pp.  1-31. 

Reference  for  additional  comparisons  of  schemes: 

R.W.  Tramel  and  R.H.  Nichols,  "A  Highly  Efficient  Numerical  Method  for 
Overset-Mesh  Moving-Body  Problems,"  AIAA-97-????. _ 
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APPENDIX  H.  GRID  MESH  OF  COMBUSTORS  USING  FAST. 


Figure  41 .  Grid  Mesh  of  2.25  inches  Diameter  Cylindrical  Combustor  (Left) 
Figure  42.  Grid  Mesh  of  2.5inches  Diameter  Cylindrical  Combustor  (Right) 


Figure  43.  Grid  Mesh  of  2.75  inches  Diameter  Cylindrical  Combustor  (Left). 
Figure  44.  Grid  Mesh  of  3  inches  Diameter  Cylindrical  Combustor  (Right). 
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Figure  45.  Grid  Mesh  of  2.25  inches  Diameter  Scalloped  Combustor  (Left). 
Figure  46.  Grid  Mesh  of  2.5inches  Diameter  Scalloped  Combustor  (Right). 


Figure  47.  Grid  Mesh  of  2.75  inches  Diameter  Scalloped  Combustor. 
Figure  48.  Grid  Mesh  of  3  inches  Diameter  Scalloped  Combustor. 
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Figure  49.  Grid  Mesh  of  2.75  inches  Diameter  Modified  Scalloped  Combustor  (Left). 
Figure  50.  Grid  Mesh  of  3  inches  Diameter  Modified  Scalloped  Combustor  (Right). 
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APPENDIX  I.  TEMPERATURE  AND  PRESSURE  CONTOUR  PLOTS  OF 

COMBUSTORS. 


Ml  IJ»  1.75  2.21  2.M  3.14  3.41  4.47  4.54  5.M 

Figure  5 1 .  Temperature  and  Pressure  Contour  Plot  of  2.25  inches  Diameter  Cylindrical  Combustor. 


n.44  n.ST  n.^  I.M  1J4  1J3  l.i2  IM  1J*D  J.lf  1.34 


nJi  D.4]  1.75  IJI  2.44  3.14  3.41  4J|r7  454  SM 


Figure  52.  Temperature  and  Pressure  Contour  Plot  of  2.5inches  Diameter  Cylindrical  Combustor. 
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TnnpEruLDrr 


UJB  I).??  I.IB  1.41  1,72  2.ni  2JS  J.H  J..M  3Ji,  3JT 


itJH2  05  1.75  02  2.«S  3.14  3.S1  4.(17  ilJ4  J.CO 


Figure  53.  Temperature  and  Pressure  Contour  Plot  of  2.75  inches  Diameter  Cylindrical  Combustor. 
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Figure  54.  Temperature  and  Pressure  Contour  Plot  of  2.25  inches  Diameter  Scalloped  Combustor. 


92 


TmperuLnrr 


1>JB  tt.6T  dJjti  l.M  1.33  IJl  l.lil  iJttt  l.Wj  2.17  2JS 


DJi  OS  1.T5  J.iJ  2.*S  3.U  3.S1  4.(17  ilJ4  J.CO 


Figure  55.  Temperature  and  Pressure  Contour  Plot  of  2.5  inehes  Diameter  Scalloped  Combustor. 


T^mpcrullurf 
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Figure  56.  Temperature  and  Pressure  Contour  Plot  of  2.75  inches  Diameter  Scalloped  Combustor. 
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Tempcrulure 


OJ  2.0  3J  5.0  6.5  8.0  9JS  11.0  12.5  14.0  15J 


0.56  0.81  UO  1.75  2J2  2.68  5.14  5.61  4.07  454  5.00 


Figure  57.  Temperature  and  Pressure  Contour  Plot  of  2.75  inches  Diameter  Modified  Scalloped  Combustor. 
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Figure  58.  Temperature  and  Pressure  Contour  Plot  of  3  inches  Diameter  Modified  Scalloped  Combustor. 
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APPENDIX  J.  TEMPERATURE  PLOTS  OF  COMBUSTORS. 


Figure  59.  Temperature  Plot  of  2.25  inches  Diameter  Cylindrical  Combustor  (Left). 
Figure  60.  Temperature  Plot  of  2.5  inches  Diameter  Cylindrical  Combustor  (Right). 


Figure  61.  Temperature  Plot  of  2.75  inches  Diameter  Cylindrical  Combustor  (Left). 
Figure  62.  Temperature  Plot  of  2.25  inches  Diameter  Scalloped  Combustor  (Right). 
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Figure  63.  Temperature  Plot  of  2.5  inches  Diameter  Scalloped  Combustor  (Left). 
Figure  64.  Temperature  Plot  of  2.75  inches  Diameter  Scalloped  Combustor  (Right). 


Figure  65.  Temperature  Plot  of  2.75  inches  Diameter  Modified  Scalloped  Combustor  (Left). 
Figure  66.  Temperature  Plot  of  3  inches  Diameter  Modified  Scalloped  Combustor  (Right). 
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APPENDIX  K.  MACH  NO.  PLOTS  OF  COMBUSTORS. 


Figure  67.  Mach  No.  Plot  of  2.25  inches  Diameter  Cylindrical  Combustor  (Left). 
Figure  68.  Mach  No.  Plot  of  2.5  inches  Diameter  Cylindrical  Combustor  (Right). 


Figure  69.  Mach  No.  Plot  of  2.75  inches  Diameter  Cylindrical  Combustor  (Left). 
Figure  70.  Mach  No.  Plot  of  3  inches  Diameter  Cylindrical  Combustor  (Right). 
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Figure  7 1 .  Mach  No.  Plot  of  2.25  inches  Diameter  Scalloped  Combustor  (Left). 
Figure  72.  Mach  No.  Plot  of  2.5  inches  Diameter  Scalloped  Combustor  (Right). 


Figure  73.  Mach  No.  Plot  of  2.75  inches  Diameter  Scalloped  Combustor  (Left). 
Figure  74.  Mach  No.  Plot  of  3  inches  Diameter  Scalloped  Combustor  (Right). 
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Mach  No. 


Figure  75.  Mach  No.  Plot  of  2.75  inches  Diameter  Modified  Scalloped  Combustor  (Left). 
Figure  76.  Mach  No.  Plot  of  3  inches  Diameter  Modified  Scalloped  Combustor  (Right). 
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APPENDIX  L.  VORTICITY  MAGNITUDE  PLOTS  OF  COMBUSTORS. 


Yorticity  Magnitude 


Yorticity  Magnitude 


I 


■ 


Figure  77.  Yorticity  Magnitude  Plot  of  2.25  inches  Diameter  Cylindrical  Combustor  (Left). 
Figure  78.  Yorticity  Magnitude  Plot  of  2.5  inches  Diameter  Cylindrical  Combustor  (Right). 
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Figure  79. 
Figure  80. 


Yorticity  Magnitude  Plot  of  2.75  inches  Diameter  Cylindrical  Combustor  (Left). 
Yorticity  Magnitude  Plot  of  2.25  inches  Diameter  Scalloped  Combustor  (Right). 
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Figure  8 1 .  Vorticity  Magnitude  Plot  of  2.5  inehes  Diameter  Scalloped  Combustor  (Left). 
Figure  82.  Vorticity  Magnitude  Plot  of  2.75  inches  Diameter  Scalloped  Combustor  (Right). 


Figure  83.  Vorticity  Magnitude  Plot  of  2.75  inches  Diameter  Modified  Scalloped  Combustor  (Left). 
Figure  84.  Vorticity  Magnitude  Plot  of  3  inches  Diameter  Modified  Scalloped  Combustor  (Right). 
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